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The physics of hquids in porous media gives rise to many interesting phenomena, including im- 
bibition where a viscous fluid displaces a less viscous one. Here we discuss the theoretical and 
experimental progress made in recent years in this field. The emphasis is on an interfacial descrip- 
tion, akin to the focus of a statistical physics approach. Coarse-grained equations of motion have 
been recently presented in the literature. These contain terms that take into account the pertinent 
features of imbibition: non-locality and the quenched noise that arises from the random environ- 
ment, fluctuations of the fluid flow and capillary forces. The theoretical progress has highlighted 
the presence of intrinsic length-scales that invalidate scale invariance often assumed to be present in 
kinetic roughening processes such as that of a two-phase boundary in liquid penetration. Another 
important fact is that the macroscopic fluid flow, the kinetic roughening properties, and the effective 
noise in the problem are all coupled. Many possible deviations from simple scaling behaviour exist, 
and we outline the experimental evidence. Finally, prospects for further work, both theoretical and 
experimental, are discussed. 
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PACS numbers; 74.60.Ge, 05.40.-a, 74.62.Dh 



I. INTRODUCTION 

It is easy to do qualitative observations on the physics 
of fluid penetration in inhomogeneous media: a drop of 
coffee on a napkin or a sugar cube held partly in the same 
coffee cup are enough to demonstrate two fundamental 
facts. A moving interface is formed between the wet and 
non-wet regions of the medium. It becomes apparent that 
the disordered pore structure and uneven surface of the 
paper or the structure of the sugar cube both influence its 
behaviour: the interface is clearly rough. Furthermore, 
the dynamics of the phenomenon slows down with time, 
meaning that the wetted area of the napkin or volume 
of coffee absorbed by the cube increases more and more 
slowly. In fact, the average position of the interface H 
of the wet front usually increases in time as H{t) ^ t^l"^ . 
The coffee drop shows an example of spontaneous imbi- 
bition, and it obeys what is known as Washburn's law 
[^39'|. The force driving the liquid from the liquid reser- 
voir to the front between the wet region and the air in 
the medium has a weaker and weaker effect on the total 
flow as the distance between these two gets larger. Such 
a naive first glance at an imbibition experiment can be 
seen in Figure ^ 

There are many similar examples of situations in which 
a liquid invades a porous medium and pushes aside an- 
other viscous liquid or gas. They are often of importance 
for technological applications or as ingredients in another 
field than the physics of fluids. In Table |l] we list some 
scenarios in which imbibition plays a role — they range 
from oil recovery (using water to displace it out of rock) 
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FIG. 1: Front of black ink sucked into a paper towel, a) 
digital photograph, ca. 1200 pixels horizontal resolution, dark 
and light grey values were enhanced to black and white, b) 
high resolution scan (1000 dots per cm) of a small part (ca. 
0.8 cm wide) in grey-scale. The structure of the medium, as 
the flbres on the paper top surface, and its effect on the fluid 
front becomes visible. 



to biology (water in living organisms) and manufactur- 
ing processes. The flow of liquids through porous media 
thus forms a very vast fleld which combines the porous 
structure of the medium with the surface chemistry and 
physics of the liquids and/or gases involved and is char- 
acterised by several parameters such as the viscosity con- 
trast of the fluids, their wettability and surface tension as 
well as the displacement rate. The simplest realisation 
of imbibition involves two immiscible fluids, one being 
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displaced by the other, both hquids being characterised 
by their viscosity 77 and the surface tension a of their 
interface. 

The sohd matrix interacts with the fluids through their 
wetting properties, as described by the Young-Duprc 
equation for a drop a hquid in contact with a sohd sur- 
face, shown schematicaUy on Fig. 12 



Cs2 — <^s\ = f cos B 



(1) 



where cr^i and Usi are the surface tensions of the sohd 
with fluids 1 and 2 respectively. The contact angle Q 
determines whether liquid 1 is wetting {Q < n/2) or non- 
wetting {6 > tt/2). 

Imbibition means that a wetting fluid displaces the 
non-wetting one, while the opposite case is called 
drainage. Spontaneous imbibition takes place when the 
invading fluid does so under the sole influence of capil- 
lary forces, with no external pressure. Forced imbibition 
involves a combination of capillary phenomena and an 
externally enforced flow rate or pressure difference. 



TABLE I: Practical and experimental realisations of imbibi- 
tion 

Oil recovery Displacement of a liquid by an- 

other, possibly in presence of a 
third phase 2, 61, 77, 241, 310, 349] 



Printing processes 

Food industry 
Biological sciences 

Surface chemistry 



Ink penetration in pa per |28al : 

Coating of paper l266l |280| : Ab- 
sorbing materials |30ll | 



Cooking [2^; Wine filtering 



Fluid transport in plants or imbibi- 
tion of water into seeds (see Section 
IIV Cll : Water penetration into soils 
|33|: Medical applications |2()^ 

Con tact ang le measurements 
|5a lil^, 113; Droplets on surfaces 



Composite materials Invasion of voids by a resin or a 
metal in filer or metal-mctal com- 
posites. UlilsO, 92, 230, 231, 232] 



Textiles 



Construction 



Behaviour of garme nts i n the pres- 
ence of liquids [H [isl Ell llS? 



Water penetration into concrete or 
cement pastes 50, 195J 



Although empirical relations for the flow of liquids 
through a porous medium existed for a long time, an ef- 
fort, in part inspired by statistical mechanics, to quanti- 
tatively understand and predict the flow led to the study 
of pattern formation or the geometry of the regions oc- 
cupied by the invading/receding fluids. At this level, the 
physics is dependent on the combination of viscous and 



capillary forces at the boundary between invaded and 
"dry" regions. The viscous part arises due to the fluid 
flow in the (partly) saturated pore space, and the relative 
importance of viscous and capillary effects is described by 
the capillary number : 



77 V 
a 



(2) 



where r\ is the viscosity of the fluid, v its average velocity 
and G the intcrfacial surface tension. 
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FIG. 2: Liquid drop in contact with a solid surface. The 
contact between the surface and liquids 1 and 2 is charac- 
terised by surface tensions Osi and (Ts2 respectively while the 
interface between the two fluids has surface tension a. If the 
contact angle Q < 7r/2, fiuid 1 wets the solid. 

The stability and existence of an interface was ex- 
perimentally investigated through the injection of liq- 
uid into various specially d esigned porous networks by 
Lenormand |187| (see also jl85l II86I ISSSf l who quali- 
tatively summarised it as a phase diagram with three 
possible outcomes, illustrated on Fig|31 (i) discontinuous 
formation of wetted domains due to surface flow in pores, 
(ii) formation of non-compact branched structures due to 
weak surface tension, and flnally (iii) compact domains 
with well defined interfaces. 

A phase diagram in terms of both the capillary number 
and the viscosity ratio of the two liquids (or liquid/gas) 
involved can then be established, illustrating the com- 
petition of viscous and capillary forces. The former sta- 
bilises the interface while the latter, if dominating, leads 
to situations like (ii) above. 

If imbibition dynamics are dominated by capillary 
forces and pore-level invasion mechanisms like film flow, 
the effective surface tension of an interface can van- 
ish. Then percolation-like phenomena can ensue, which 
means that the medium is apparently penetrated in dis- 
joint clusters of the imbibing liquid. On the pore scale, 
one thus has to deal with either piston-like displacement 
or such gradual processes, leading flnally to "snap-off" , 
as narrow pore throats b ecome completely filled by the 
invading liquid p39ll283|. 
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FIG. 3: Various cases of a fluid (oil) displacing air in a net- 
work, with varying pore size distributions. From left to right: 
a square network geometry, a triangular one, and finally a 
square one with very wide pore sizes. As the Capillary num- 
ber Ca is varied, the effective surface tension of the cluster 
of invaded pores changes. Note the results for logCa = — 6 in 
particular (after Lenormand, 187']). 
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FIG. 4: A phase diagra m for the geometry of imbibition as 
proposed by Lenormand }l8j| . with the capillary number Ca 
and the viscosity v being varied on one hand, and the pore- 
level geometry on the other hand. Notice the possibilities of 
apparently non-local (discontinuous) invasion, and of vanish- 
ing surface tension. 



If the imbibing fluid is sufficiently more viscous than 
the defending one the invasion front is compact. Typi- 
cally, as in the example of the sugar cube or the paper 
napkin, it is not flat but presents a random rough struc- 
ture. This additionally complicates any effective theory 
or equation of motion for the invading fluid. 

The main issues to be discussed in this review arti- 
cle are related to the rough interface between the two 



phases involved in a compact imbibition invasion. This 
is accomplished by coarse- graining from the level of in- 
dividual pores to the average position of the interface, 
meaning that the micro-structure of the medium is aver- 
aged out when possible. This is particularly complex in 
the case of spontaneous imbibition and can fail, due to 
the slowing down of the fluid penetration and/or due to 
the presence of competing mechanisms in the pore scale 
invasion dynamics. 

This focus is also of fundamental interest to any im- 
bibition process since a correct interfacial description re- 
quires the knowledge of the time- and length-scales that 
control the entire process. For many practical applica- 
tions like soil mechanics or oil recovery it is of importance 
to understand this kind of "upscaling" ||29||, l24lL l336| | in 
order to estimate remaining non- wetting fluid saturations 
or relative permeabilities from laboratory scale measure- 
ments or from microscopic simulations. Future advances 
can also be expected in more complicated scenarios, as 
in the case of non-Newtonian fluids, in the presence of 
chemical reactions between the fluids involved, as well as 
in the rapidly developing field of microfluidics. 

The process of imbibition is affected by noise, part 
of which stems from ever-present thermal fluctuations, 
while the rest is due to the quenched, frozen-in struc- 
ture of the porous medium. The quenched nature of the 
noise becomes particularly important if the phase inter- 
face moves slowly in avalanche-like behaviour. Then the 
dynamics consists of localised bursts whose description 
directly couples the noise and the dynamical fluctuations 
of the interface. 

An interesting theoretical question is then the univer- 
sality of these phenomena, i.e., whether the statistical 
description of the interface is similar to the roughness 
observed in flre fronts, cracks and rupture lines, domain 
walls in fcrromagnets. The same question arise with 
respect to the description of burst and avalanches in 
connection to the concept of Self-Organised Criticality. 
These questions are far from academic. In one part, the 
description of roughening in non-equilibrium phenomena 
has received lots of attention from the theoretical side 
(witness the large number of publications devoted to the 
Kardar-Parisi-Zhang (KPZ) equation [14111). but it is only 
recently that experimentalists were able to convincingly 
demonstrate the link between the propos ed th eoretical 
models and the observed phenomena (see j229 ] and ref- 
erences therein for KPZ behaviour in fire fronts). The 
same is true for roughening in imbibition. 

Although it was examined in the early 1990's '?,'45'| in 
connection with percolation theory and deviations from 
KPZ behaviour, recent work [s^ l99|| pointing to the im- 
portance of fluid conservation, has led to a flurry of new 
experimental results 101, 122, 311, 312, 313J supported 
by further theoretical work jl78Ll26lj . Even though many 
details remain obscure, we feel that the general theoreti- 
cal picture of roughening in imbibition is now well estab- 
lished. The goal of this review is thus to highlight, from 
a statistical physics point of view, the central aspects of 
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imbibition that are understood, to relate them to exist- 
ing experimental results and to point to remaining gray 
zones that deserve further studies. At the same time, we 
want to contrast this with the array of experimental evi- 
dence and numerical simulation models used for applied 
purposes and in connection to multi-phase flows. New 
theoretical results on the dynamics of avalanches and on 
interface roughening at constant flow and for columnar 
disorder ar e also presented an d related to the recent ex- 
periment [HI llS 113 in Sections imniinini 

and lIIIG6l 

The experimental situation is reviewed after that, with 
a brief remainder of possible complicating effects and 
comments on non-Washburn-like scaling in various se- 
tups. Next follows an account of various experiments on 
front roughening, ranging from Hele-Shaw cells to pa- 
per imbibition, divided between the forced fluid flow and 
the spontaneous cases. In contrast to the recently fo- 
cussed theoretical progress, the experimental picture is 
much more delicate and we discuss it in the light of the 
theoretical sections. It turns out so that there are many 
confusing results, some of which can be related to the 
theoretical ones, and some of which clearly call for fur- 
ther work. 

It is now clear that the physics of imbibition in dis- 
ordered media couples randomness, kinetic roughening 
properties^ and interface velocity in an complex fash- 
ion. This is illustrated schematically in Figure |5| The 
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FIG. 5: Schematic representation of elementary processes in- 
volved in quasi-two-dimensional imbibition. 

randomness acts as an effective noise that roughens the 
interface, but the same randomness also influences the 
fluid flow in the bulk, and thus the interface velocity and 
the net volume flow in the imbibition process. There- 
fore, we must conclude that all these aspects are inti- 
mately coupled and cannot be discussed separately. This 
creates intriguing experimental and theoretical questions 
and brings engineering and technological interests, sur- 
prisingly, close to those of fundamental statistical me- 
chanics. 



At the basis of models and experiments dealing with 
imbibition interface roughening lies the understanding 
of bulk fluid transport and saturation in the invaded 
medium. In the section to follow next, an overview is 
given from basic hydrodynamics in the medium to spe- 
cial effects related to non-Newtonian fluids. 



II. PHYSICS OF IMBIBITION 
A. Fluid flow in porous media 

The detailed description of liquid flow through a 
porous medium is greatly complicated by the many time 
and length scales that are involved. The flow takes place 
at the level of the pores (say 1 /im or less) but the quan- 
tity of interest is the flow averaged over the whole macro- 
scopic sample. 

Similarly, a pore may be invaded vary rapidly or may 
remain blocked for the whole invasion time, which may 
be several hours or even days in spontaneous imbibition. 
A description of the flow makes sense only over a "rep- 
resentative volume clement" , loosely defined as the min- 
imal volume that can be defined such that the properties 
of the flow (and of the porous medium itself) within it 
remain statistically similar no matter w here it is placed 
in the bulk of the medium (see |227l l348| | for some recent 
advances using tomographic techniques). 

At this level, it was already empirically found long ago 
by Darcy that the flow is essentially proportional to the 
pressure gradient across the medium, a result which can 
intuitively be explained by considering the flow of a liq- 
uid through a capillary tube. A consequence of this result 
is that the height iJ of a fluid column invading sponta- 
neously a porous medium from a reservoir increases in 
time as H{t) ~ t^/^, due to a combination of fluid con- 
servation law (since any amount of invading fluid must 
be transported from the reservoir) and capillary forces. 
Before going further into the concepts of interface rough- 
ening, it is interesting to see how this result arises, and 
under which conditions it can be expected to be valid. 



1. The Lucas-Washburn description 

The simplest way to illustrate imbibition is the capil- 
lary rise. It translates easily to the basic phenomenolog- 
ical description of flow in porous media and it represents 
an important microscopic mode of flow propagation. As 
illustrated in Fig. El we consider a capillary tube of length 
L and radius i?0 7 immersed into a reservoir at ambient 
atmospheric pressure Pq. 

The two fluids are immiscible with an interface, of sur- 
face tension a, located at a height z = H above the reser- 
voir and we assume that the fluid in the region z < H 
wets the solid walls. A meniscus is formed at the inter- 
face, characterised at equilibrium by the contact angle 9 
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FIG. 6: Illustration of capillary rise, the fluid at z = is 
in contact with a reservoir at atmospheric pressure Po and 
the top of the tube is either closed or open. The fluid from 
the reservoir wets the walls of the tube, which means that 
a curved interface, to which corresponds a capillary pressure 
difference Pc exists across the interface. 



obtained from Eq. . The combination of the curved in- 
terface with the surface tension creates a pressure differ- 
ence or capillary pressure pc = AP = 2a cos 9 /Rq. This 
Gibbs-Thomson or Laplace effect gives rise to a thermo- 
dynamical force that will move the fluids. It also applies 
in more general imbibition problems directly at the level 
of individual pores. 

In its full generality the solution of the flow dynamics 
for this problem needs to take into account the detailed 
flow right at t he contact line between the fluids and the 
solid |70l . Il43| as well as inertial effects at the entrance 
of the tube |319l | . For now, these effects are neglected, 
but are discussed in greater detail below. Neglecting the 
structure of the meniscus means that the pressure field 
changes only along the length of the tube, P = P{z), 
with the interface described mathematically by the height 
H. Likewise, neglecting inertia means that the fluids are 
only descr ibed in terms of their viscosity rj using Stokes' 
Equation ^7^. 



i]\7'^v = VP 



(3) 



Three different cases can be then be easily examined: 

(i) Two incompressible Fluids: The two fluids are de- 
scribed by Stokes' equation, with a longitudinal velocity 
field depending only on the radial coordinate v — z Vz{r): 



1 d 
r dr 



dv,,{r)\ dP{z) 



dr 



dz 



(4) 



where i = w or n denotes respectively the wetting (z < H) 
and non-wetting (z > H) fluids. The tube is open at 



the top, to allow the second fluid to leak out. The 
incompressibility condition leads to a Laplace equation 
V^P — for the pressure field, where the instantaneous 
position of the interface appears only as a boundary con- 
dition 



P{z = 0) = Po 
P{H+)-P{H-)=p, 
P{z = L)=Po. 



(5) 



The Laplace equation for the pressure means that the 
pressure gradient is constant, so that each fluid has a 
Poiseuille velocity field 



1 

4?7j 



(6) 



and the velocity of the interface is associated with the 
average flow velocity 



ttRI 



dH 
It 



da.Vz^{r) = / da.Vzn{r) 



(7) 



where da is the area element of the tube, which leads to 
a pressure 



P{z<H) = - 



P{z>H) 



Vn 



T]^H + T]n{L~H' 



Vn 



Pc Z + Po 



(8) 



and to an interfacial velocity 

Pc 



Pc{z ~L)+Po 



dH 
~dt 



Rl 



8 T]^H + T]niL- H) 



(9) 



The wetting fluids thus continuously displaces the non- 
wetting fluid until it occupies the whole length of the 
tube. Notice that at the beginning, it does so as dH/dt ~ 
Pc/H, which implies an initial motion of the interface 
H{t) ^ 

(u) Second fluid is compressible: If the second fluid is a 
very compressible liquid or a gas, it simply adjusts itself 
to the pressure of the wetting fluid and does not sup- 
port any pressure gradient. In the case of an open tube, 
the second fluid is the ambient gas, and its pressure is 
constantly at the atmospheric pressure Pq, such that 



P{z<H) 



-P.-+PO 



(10) 



P{z>H) = Po 



In this case , the classic result of Washburn |339j and 
Rideal |278l | is obtained: the interface position is de- 
scribed in time as 



RIpc 



[t - to) 



(11) 
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where H{to) is the initial height. This result is also ob- 
tained for two incompressible fluids (Eq. as long as 
r]njH > r]nL. 

If the tube is closed at the top, so that the gas cannot 
escape, the pressure of the gas phase increases as the fluid 
occupies more and more volume. If the gas was initially 
filling the whole tube with pressure Pq, then the pressure 
if the interface is at height H, 



P{z >H)=Po 



(12) 



where 7 — Cp/Cy, the ratio of the heat capacities at 
constant pressure and volume respectively. The flow field 
of the wetting fluid is not modified, and the resulting 
interface equation 



dH_ ^ Rl 1 P{H) -Po-Pc 
dt 



8 77t, 



H 



(13) 



Again, as long as P{H) is not too far from the atmo- 
spheric pressure, the Washburn-Rideal result is obtained. 
However, in real porous media, trapping of gas by the in- 
vading fluid may occur, in which case the propagation of 
the menisci, given by Eq. I|13|) may be markedly different 
from Washburn behaviour Eq. (|11|1 . 

(iii) Capillary rise and gravity: gravity acts on the flu- 
ids through their density p as 



dP{z) 
dz 



Pi 9 



(14) 



where g is the gravitational constant. Assuming that the 
second fluid is the ambient gas, with an open tube, the 
interfacial rise is 



dH_ 
~dt 



K 



Hen 



(15) 



where we define the permeability k = and the 

equilibrium height H^q = Pc/{p9), at which the hydro- 
static and capillary pressures are balanced. The relevant 
parameter is the equilibration time Tcq = Hcqri/Kg{p)^. 
For times t <C Tgq (which corresponds to H ^ Hcq) or 
in absence of gravity ( "horizontal imbibition" ) the time- 
dependence follows the square-root law, Eq. ((TT|l seen 
earlier At later times, t ^ Teq, the finite equilibrium 
height Hcq (if it exists) is approached exponentially 



H{t) ^ H,q (1 



(16) 



in a way that is simply related to the equilibrium quanti- 
ties Hcq and Tcq. This simple law is interesting in many 
aspects since the relevant parameters (contact angle, sur- 
face tension) are easy to determine and the resulting 
easy-to-grasp predictions make it attractive and suitable 
in interpreting experiments. 

However, these results may be invalid for a number of 
reason s. Fo r example, in a liquid-liquid system Mumley 
et al. [242j |. comparing dry and prewetted tubes, have 



shown that viscous dissipation at the contact line and 
precursor film dynamics 7Q . il4S | may lead to both Wash- 
burnian or slower (in terms of the rise exponent) be- 
haviour depending on the contact angle and the condition 
of the capillary rise. 

Inertial effects as Bosanquet-flow [s^ may also be im- 
portant in the early stage of pore inva sion, before the 
dynamics is described by Eq. Hll|l |319l |. Dimensionally, 
it is easy to see that the characteristic time scale over 
which inertial effects will be important ~ pR^/r], after 
which the usual Washburn dynamics follows. Before this 
time the fluid enters the capillary as 



H{t) 



— 

\pRo 



e/2 



t. 



(17) 



— 1/2 

Although the time n can be very small, the Rq de- 
pendence of the rise on the radius of the capillary must 
be compared to the Rq of Eq. l(T^ . Although the ef- 
fect of inertia lasts for a very short time in small pores, 
the capillary rise can nevertheless be quite rapid. Exper- 
imentally inertial effects can give rise to effective initial 
height and time, Hq and to in Eq. H15|l that must be taken 
into account in fitting measured data |l73l Il74| . 



2. Macroscopic description of flow in porous media 

Equation H15|l is also often used in more general terms 
for porous media. Already at the level of Stokes' equa- 
tion, a simple dimensional analysis implies that v ~ 
(k/?7)VP where k must have units of (length)^. In fact, 
a common description of si ngle- phase fluid flow in porous 
media is Darcy's equation [222, [2£ 



--(VP-pg), 



(18) 



where (q) is the average volume of fluid transported per 
unit time per unit cross-section of the porous medium, 
and K and 77 are the average permeability and viscosity re- 
spectively. This equation is obtained by coarse-graining 
the structure of the porous medium, ignoring all fluctua- 
tion and pore scale effects, using a line of reasoning based 
on "representative volume elements" Q, 295]. 

The extension to two-phase or multi-phase flow is in 
principle straightforward, the dimensionless saturations 
of wetting and non-wetting fluids Sw and s„ are deflned 
as the ratio of fluid present with respect to void space 
in a given representative volume element. Alternatively, 
the concentrations are defined as Ci = Vsi where V is 
the porosity. The dynamics of the concentration is then 
obtained through continuity equations to refiect the con- 
servation of liquid. If the fiuid densities and the porosity 
are constant, 



r 

V 



~dt 

dSn 

dt 



(19) 
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together with the constraint s^, + s„ = 1. The flux 
of each fluid obeys Eq. (|18|) and a coarse-grained capil- 
lary pressure is included as a closure relation between the 
pressure in the two fluids Pw ~ Pn — Pc- Since the per- 
meability is only a geometric factor, both fluids should 
have the identical values of k and the pressure gradient, 
VP ~ z pc/H{t) then follows immediately. 

If the imbibition front remains compact, the volume- 
flow per unit area is nothing but the velocity of the inter- 
face, which leads to (q) = dH/dt and to the results de- 
scribed in Eq. lO, and Eq. if the second fluid cannot 
sustain a pressure gradient. This however assumes that 
the fluid pushed away has an easy escape, a situation 
denoted in the petrophysics community as "co-current 
flow" . 

What kind of dynamical behaviour can now be ex- 
pected in realistic porous media? Simplistic analysis pre- 
dicted by the above Washburn law, Eq. can be inval- 
idated by many factors. For example, evaporation (to be 
examined in more details in Section flll Gl) may slow down 
imbibition when the quantity or volume of liquid needed 
to follow the typical Washburn scaling is not available. 
Similar effects arise if the solid phase matrix is perme- 
able and acts as a sink. Likewise, for three-phase flows, 
often important in industrial problems, the separation of 
the two liquids may cause the introduction of new time- 
scales. In general, simply having several length-scales to 
describe the structure of the porous medium can lead to 
complications that obfuscate the Washburn scaling. In 
the following, we briefly discuss when this scaling may 
fail, or what could be expected in various kind of porous 
media. 

Despite all these possible complications Washburn be- 
haviour is often observed experimentally. It is therefore 
an useful approach to study roughening first under well 
controlled situations that lead to the t^^^ progression of 
the average interface, and then to look at possible de- 
viations. We will therefore devote the major part of 
this Review to Washburnian imbibition. The expected 
roughening behaviour when H{t) ^ with (5 7^ 1/2 is 
commented on at special occasions, in particular in the 
concluding Section. 



3. Pore geometry and inertia 

At the microscopic level, it is clear that any simple 
picture of a porous medium as a collection of capillary 
tubes with effective permeability k, is wrong, the pore 
geometry can be extremely important, in particular when 
sharp corners or edges are present. 

It is already established that a compact (albeit rough) 
imbibition front is formed if the viscosity of the invad- 
ing fluid is sufficiently greater than the viscosity of the 
defending fluid, as indicated on the phase diagram pro- 
posed by Lenormand 187], reproduced in Fig.^ Thus, 
although containing interesting physics at the level of 
pattern formation, the question of disconnected inva- 



sion clusters is trivially not relevant in the context of 
interface roughening (see e.g. refs. 0; Il26j for perco- 
lation approaches). The main exception to this rule is 
the possibility of a "dynamic" transition between well- 
defined interfaces and fingering instabilities in sponta- 
neous imb ibition due to the slowing down of the interface 
|l2Cl l328l | . A second possibility is the role of prewetting 
layers, whose presence may again ren der an interfacial 
description meaningless |20ll l202l |203| | , since the essen- 
tial dynamics are ruled by the prewetting. 

In most cases it cannot be decided a priori whether 
the expected macroscopic Washburn-Darcy behaviour is 
valid or not. The early studies of Lenormand, later aug- 
mented by by the studies of Bernadincr, Knackstedt et al. 
and others, have managed to shed further light into the 
question of relevance or irrelevance of pore-scale physics 
I25I l302j| . Three basic invasion processes were observed. 
A cylindrical pore may be invaded either as in the basic 
capillary rise (piston fiow) or by surface flow, followed by 
collapse (snap-off) of the film, as shown in Fig. |7| At 
the crossing of several cylindrical pores, there is a dras- 
tic change in the shape of the meniscus, followed by a 
rapid pinch-off leading to the invasion of the pore (Fig.|SJ). 
Each of these modes of invasion are characterised by dif- 
ferent length scales, which can lead to a dynamical be- 
haviour that is very different from the one predicted by 
Eq. Hll|) . Blunt et al. [s^ argue that the width of a film 
fiow zone in networks (under forced imbibition) should 
scale as ~ {wt/Ro) where Wt/Ro is a rough esti- 
mate of the ratio of tube wall roughness to its diameter. 
Percolation-style arguments for snap-off region sizes indi- 
cate that the associated length-scale increases much more 
weakly with 1/Ca 113 . 

One question related to the pore-level description is 
whether correlations matter. In various types of rocks 
the local permeability may be correlated over consider- 
able length-scales due to geological processes. The ev- 
idence from numerical modelling based on microscop- 
ically faithful descriptions of pore-networks seems to 
indicate, that multi-phase flow phenomena as imbibi- 
tion in particular are much more prone for correlation- 
induced effects than, e.g. , simp le on e-phase p ermea bility 

Elii, EE il im EHi ill I253 iHl Hoi HH ■ 

Typical quantities where this would be visible are rela- 
tive permeabilities for given saturations, and remaining 
saturations of the non-wetting fluid. Note that there are 
no studies of the front dynamics in imbibition in the pres- 
ence of "typical" correlations of pore structure, as arising 
in empirical contexts. 

In addition, inertial effects |28lll299j may be important 
in the early stage of pore invasion, where they can lead 
to a preferential invasion of small pores (see Eq. (|17|) ') 
281, 299]. But inertia can also be important at all stages 
of the capillary invasion, whenever a sudden change in the 
fiow occurs, as in film snap-off, or pinch-off at a corner. 
In thes e cases, th e front position could advance as t'^ with 
k < 1 |27Ct l27lj . On short time scales (smaller than a 
few seconds), the front penetration may also proceed in 
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FIG. 7: In large pores or channels, the fluid can fi rst c over 
the surfaces and then collapse to flU the whole void |l87j| . 





a) 



b) 




c) 



d) 



FIG. 8: Examples of pore invasion mechanisms |l87j : the 
progress is gradual from a) to b) as the meniscus radius in- 
creases, and leading, at c), to a situation when the fluid sud- 
denly invades the channels. 



pores first with the establishment of a prewetting layer 
through diffusion 90., 114] or (precursor) film How 



4-. Dynamical Saturation 

In general capillary flow occurs simultaneously with 
film flow and can lead to a gradient in the concentration 
of the invading fluid: complete saturation is established 
only over a certain distance. In addition, trapping of the 
defending fluid in certain regions modifies the flow of the 
imbibing fluid and may block some paths of progression, 
leading to a reduced permeability for the invading fluid 
(e.g., Eq. fl^ for the simplest case). It is also often 
possible that the imbibing liquid interacts directly with 
the porous matrix, changing directly the spatial structure 
of the porous medium and leading to a slow change i n th e 
permeability (as observed in concrete with water lli^ 
or in ordinary paper [s^l. 



These phenomena can, to some extent be taken into 
account by introducing different permeabilities for each 
fluid in Eq. ifH^ : = Ki(VPi — pigz)/rii with i = w,n. 
In addition, the permeabilities and capillary pressure will 
often be function of the saturation of the medium, Ki = 
Kj(sj) and pc = Pcisw)- 

However, even this basic set of equations has to be 
augmented to account properly for surface tension and 
surface areas |l23i Il24j . In fact, the two- fluid interfa- 
cial area is one of the fundamental quantities — it also 
partly defines the driving thermodynamic force in the 
case of spontaneous imbibition. Many authours have 
addressed this question by network models and exper- 
iment s on one hand or by th ermodynamical considera- 
tions [lOd llOSl I255I l276l |294| . Again the physics at pore 
level is complicated by the presence of precursor films 
and the question of how entrapped volumes of the non- 
wetting fluid behave. 

If both permeabilities are equal and independent of 
saturations, and likewise for the capillary pressure, then 
the ideal case of two incompressible fluids in a capillary, 
Eq. (jnj, is recovered. As mentioned above, this implies 
that the non- wetting fluid can easily be pushed away by 
the wetting liquid. This corresponds to "co-current flow" 
and — q„. The opposite case is "counter-current 
flow" , the mass-flow of the expelled liquid takes place 
in the exact ly op posite direction from the invading one: 
QiL. = — |l96j . In both cases, if the geometric setup 
allows for a one-dimensional equation we arrive at the 
effective diffusion equation 



dt 



dx \ ^ dx 



(20) 



where Mu,{krw, km, dpc/dsw) is an effective mobility and 
the subscript r refers to relative permeabilities. This has 
naturally a scaling solution with x/^/t as a scaling pa- 
rameter. Notice that this is also related to the piston-like 
fronts in the Buckley-Leverett theory 113 1 ^^"^ that 
the implied scali ng properties are of practical interest, 
as well |206l l347| . If the non- wetting fluid cannot sup- 
port a pressure gradient (and we assume that it is free 
to escape), then P„ — Pq and Richard's Equation |277j 
appears as a special case of Eq. (1^ . 

When dynamical saturation is present, the concept of 
an interface can easily become ill-deflned if there is no 
sharp jump in concentration between the wet and dry 
parts of the medium. In many cases the permeability 
is a rapidly increasing function of the saturation (e.g., 
n{sw) ~ exp(/3(su; — So)) (100 . 236J) so that a sharp di- 
vision exists between the invaded and non-invaded part 
of the medium, which also leads to the basic scaling 
H{t) ~ t^^^ for spontaneous imbibition. 

Hysteresis is also an omnipresent feature |^ . Partly it 
occurs between what commonly is called "primary" and 
"secondary" imbibition. Such history effects arise if there 
is already a presence (residual saturation) of the wetting 
fluid. 
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Saturation can be studied dynamically since the grad- 
ual increase in the local fluid volume fraction can be fol- 
lowed by Nuclear Magnetic Resonance (NMR) and by 
X-ray Computer Tomography (CT) techniques. These 
have reached already accuracy levels in which individual 
pores can be monitored as filled or not, and perhaps in 
the near future the interfacial features will become acces- 
sible as well (see ,341.] for a hydrology example) to allow 
studies of interface dynamics in detail j28, 320]). One 
of the important practical questions is the dynamics of 
entrapped fluid. Though percolation-based descriptions 
of su ch residual clusters have existed for quite some time 
|342j| . the physics of such "ganghons" is still quite open, 
often accessible only via numerical models. 

Imaging techniques are sufficient to demonstrate phe 



nomena like swelling of the solid volume fraction 32C 
in the case of disordered fibre networks (see also . 
IT3I I251I l337j 'l. Swelfing, coupled to the fact that the 
porosity undergoes a simultaneous decrease, leads to re- 
duced fluid flow towards the interface and possibly to 
deviations from Washburn behaviour. The problem of 
deformable porous media is difficult to understand in gen- 
eral terms, and involves from the theoretical viewpoint 
the solution of coupled elastic and fluid mechanics prob- 
lems (see e.g. 268, 309]). Other NMR studies have been 
recently performed in situations that are of relevance for 
constructio n en gineering (water imbibitio n into cem ent 
pastes 113, Ellj) and the oil industry 0, Ell |34|]. In 
the latter case, the displacement of oil by water is an 
important practical question, as is the achieved level of 
saturation. NMR imaging has now developed to a level, 
where one can consider the saturation vs. pore sizes — 
and distinguish between co- and counter-flow due to the 
importance of fllm flow in the latter one 55]. 

By simply comparing the mass intake and the optical 
appearance of a paper sample, it has been noted that the 
saturation and the actual visual front may both display 
similar temporal scalings (Washburn-like). This would 
define a widening "saturation front" as the partially sat- 
urated volume behind the interface that drives the imbi- 
bition process (27j . The principal question is whether the 
saturation takes place in the wake of a "front" or whether 
the process as a whole obeys diffusion into, e.g., small 
pores first |ll4j . The first limit corresponds to piston- 
like advance (a term using a one- dimensional analogy), 
and the second to a flat saturation profile (which changes 
by the time-scale contained in the scaling variable). 

Recent imaging data 0, l227l l228j show convinc- 
ingly the capabilities of both CT and NMR approaches 
in demonstrating both limits, and the cross-over in- 
between. Figure El exhibits a clear exampl e of a n exper- 
imental front with only partial saturation |349| . In the 
close future one would expect that su ch te chniques yield 
much more empirical evidence [t^, 113, l296j by looking at 
saturation profiles and concomitantly capillary pressure 
curves as the saturation is varied. 

The transport of liquid into the less-saturated regimes 
can further be complicated by the presence of an initial 




FIG. 9 : Sa turation profile of the liquid in a fibre assembly 
(plug) I320II , as a function of time and distance from the liquid 



reservoir. 



(residual) saturation |l4lj . Figure [T51 shows an example 
of how the presence of w ettin g fiuids is supposed to influ- 
ence the front geometry '136^ . The width of the front can 
spread in time and exhibit what is called "hyperdiffusiv- 
ity" via the dependence of the k{sw) on the local wetting 
fluid saturation. This is in fact related to the question 
of diffusive hydrodynamic spreading of tracers, and has 
been debated in the literature §7i, |6S] . Recent imag- 
ing experiments indicate that the hype r-dispersio n can 
be related in a power-law fashion to Sw j227l l22^ . This 
would result from saturation dynamics following film be- 
haviour in the pores. 

Further such complications become evident if the de- 
scription by Richard's equation is modifled by, e.g., a 
time-depe ndent porosity jl97i) . or if anomalous diffusion 
is invoked |l6£^ . both resulting presumably from liquid- 
porous matrix interactions. The usual saturation dynam- 
ics, contained in simulation models 283, 324] , would indi- 
cate diffusive intake — which then would exhibit similar 
scaling as the spontaneous imbibition front position |l68l ] . 
The saturation dynamics can thus omplicate the dynam- 
ics at the front due to the conser vatio n of the liquid that 
a sample intakes per unit time |305| . an issue that has 
not been studied systematically to our knowledge in spite 
of some "snapshots" in the literature, as Figure [T5I 



5. Surfactants, additives and other non-Newtonian effects 

Finally, additives to a carrier fluid, such as dye, give 
rise to phenomena that can cause a deviation from 
Washburn-like behaviour. The dye particles diffuse as 
usual Brownian particles but are also carried by the im- 
bibing fluid, which gives rise to hydrodynamic dispersion 
|292j . The dynamics of dispersion is of fundamental in- 
terest in the imbibition context for two main reasons: 
In some cases, the dye front is (maybe somewhat dan- 
gerously) identified with the "elastic interface" for which 
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FIG. 10: The swelli ng of fibres, as a function of time and 
distance, in a plug (' |32C| 'I. A concomitant change in local 
permeability, and thus fluid flow, is expected. 
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FIG. 11: The front and the mass (of the silicone fluid) in a pa- 
per sample vs. time (Washburn-scaling, ~ t^^^), according to 
Ref . . The implication is that the saturation is incomplete 
behind the front. 



theoretical models are built. The dynamical contact an- 
gle at the advancing front is generally dependent on the 
local dye concentration. This may act as a surfactant, 
with a time-dependent concentration. Similar behaviour 
can be observed for various mixtures of liquids: time- 
dependent changes to the composition vary the viscosity 
and affect the time-scales [T^I- Alternatively the absorp- 
tion can be slow on the time-scale of the dynamics 325]. 
We do not want to deal extensively with the effects that 
these complications may have on the reformulation of 
the interfacial theories for imbibition, but list below only 
some possible effects that one should keep in mind. 
If the surfactant transport is diffusive, one gets ac- 




FIG. 12: X-ray computer tomography of spontaneous 
coun ter-current imbibition of water into n-decane in diatomite 
|349| . Notice the slowing-down of the sharp interface between 
dark and light (the latter denoting water-saturated regions). 





(b) 



FIG. 13: Dynamic simulations of imbibition, showing how a 
front can exist even in conditions of "s econdary imbibition" or 
initial wetting phase saturation. |l3fil |. The capillary number 
is Ca = 10~^, and the wetting saturation 0.4. In the top 
panel, a snapshot, while the two lower ones show the relative 
permeabilities for two and three dimensional systems (128^ 
and 20'^, the contact angle 9 — 50 degrees. Initial wetting 
phase saturation was chosen as 0.063. 



cidentally Washburn-like dynamics but anomalous dif- 
fusion can clearly change this and lead to almost arbi- 
trary temporal scalings. At the very least, the expected 
length a nd time scales will change, sometimes drastically 
|56Lll27| . The effec t of t he surfactant (or time-dependent 
composition) (e.g. |205l i273ij ) can be understood based 
on, e.g., viscous dissipation at the contact line, which 
changes the dynamical contact angle. This is a research 
field in its own, related to the physics o f thin films and 
droplet spreading (see [68Ll217ll315ll316| among others). 

Generally, deviations from Washburn behaviour can 
easily exist. This is particularly s o at early times, due 
to the dissipation-induced changes |270| . The fast move- 
ment of menisci, coming from microscopic geometric con- 
siderations, coupled to the viscosity of the fluids involved, 
leads to dissipation at the front (viscous stresses) and 
should persist in all regimes and at all times in sponta- 
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FIG. 14: The rise of a column of water (and the effective dy- 
namic contact angle) vs. time in a capillary rise example. No- 
tice the changing angle and the deviations from Washburn's 
law (best fit with the dashed line) |118| . 



neous imbibition |328l l329l | . Some progress in accommo- 
dating dynamical phenomena can be made by modify- 
ing directl y the contact angle in Washburn-like effective 
equations |ll7j . see Fig.lCT 

It should not be forgotten that the porosity can be also 
be changed in a time-dependent way by the deposition of 
any particles carried by the invading liquid ^97, 98, 340], 
which again may lead to some erratic behaviour for the 
front. We note that the same conclusions can coinci- 
den tally be o btained in the case of non-Newtonian fluids 
[60ll262l|343j |. Only recently there have been some exper- 
iment al de velopments on the scale of network model sys- 
tems |327| . Similarl to three-phase systems [z^ such 
work highlights the upscaling problems when the pore- 
level behaviour varies widely. For shear-thinning fluids 
the local fluid response inside a single pore depends on 
the flow rate, i.e., on both the global front dynamics (due 
to fluid conservation) and on pore details such as size or 
connectivity. A possibly general trend is a less distinct 
front due to microscopic flows. Similarly, it should be 
noticed that ink, often used in imbibition experiments 
behaves rheologically as a non-Newtonian fluid 0| 

The treatment of microscale dynamics, such as inertia 
and film flow at the pore-level, is very delicate at the 
level of h ydrodynamic equations and initial conditions 
|l59lll93l |. However, the central point is that these details 
enter the hydrodynamical description in a coarse-grained 
way. This leads to a macroscopic front velocity which has 
either a Washburn-like behaviour, but with effective per- 
meability and capillary pressure, or to a completely dif- 
ferent scaling behaviour |302| . Microscopic simulations of 
capillary imbibiti on im ply that modified Washburn-type 
equations ensue |219j . Nevertheless, on a macroscopic 
scale, averaged over a history of penetration and a rep- 
resentative volume, the situation is of course open. 



B. Non-equilibrium processes and roughening of 
fronts 

The front between the wet and dry regions roughens as 
it propagates. This is the main emphasis in our review. 
The spatial and temporal statistics of the interfaces is 
easily observed by the naked eye, as in Fig. ^ and they 
remind, in general terms, of other "rough" objects such as 
fractures in inhomogeneous media, fronts in slow combus- 
tion, flux lines in disordere d superconduc tors, and Brow- 
nian paths in diffusion |2l[lll IHl IHIl . The existence 
of so many analogies poses the fundamental question of 
universality. If the roughness is scale- invariant^ it can be 
described by fluctuations and correlations that possess 
power-law behaviour, reflected in typical quantities such 
as the spatial and temporal two-point correlation func- 
tions, 6*2 (r) and C2(t) respectively (see Eqs. (|17j) and 
(|49|l ). These functions exhibit dynamical scaling that es- 
tablishes critical exponents, understood in the usual sense 
of statistical mechanics. The exponents — like tho se of 
the celebrated Kardar-Parisi-Zhang (KPZ) equation |l47j 
— are related to a particular universality class with its 
symmetries embedded in a Langevin-like equation for the 
interface. In Sections UTTTTTl to ImTTSl we review the ap- 
propriate parts of the theory of kinetic roughening. 

In the late 1980's much effort was put into the physics 
of fluid fronts (e.g., in disordered Hele-Shaw cells) to pro- 
duce clear power-laws with well-established exponents re- 
lated to a particular universality class. The hope in the 
imbibition context was that the non-local nature of the 
problem could be neglected in favour of a local descrip- 
tion, partially for simplicity but also because the compli- 
cations involved would otherwise seem unsurmountable 
|l57j . In reality, the experimental picture turned out to 
be much more complicated and often not directly related 
to any of the expected universality classes and the ques- 
tion of an effective theory for imbibition front remains. 

The main problems that a theoretical description must 
solve are the nature of the noise acting on the interface 
and the effect of non-locality induced by fluid conser- 
vation. At the simplest level the fluid pushing the inter- 
face flows according to Darcy's law, which already implies 
that the advances of sections of the fluid front are cou- 
pled. Often a length scale along the interface emerges, 
which is related to the average motion of the interface and 
sets the maximum extent of "local fluctuations" . 

This is intuitively clear for spontaneous imbibition 
since fluctuations have a larger instantaneous velocity 
when behind the front compared to when that are ahead 
of it. Consider an effective surface tension 7*, repre- 
senting the macroscopic energy cost of a curved inter- 
face. The curvature arising from fluctuations changes 
the pressure by AP = ^*W/£,l. (the Laplace pressure or 
Gibbs-Thomson effect) where W is the vertical and 
the lateral scale of the fluctuation. This pressure change 
slows down advanced parts of the interface. Comparing 
it to the overall pressure drop across the same vertical 
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distance, AP = pcW/H, yields a lateral length scale 

^ ^H/p,. (21) 

This indicates that fluctuations on a scale larger than 
are suppressed by the pressure field gradient Pc/H . For 
an interface slowing down in spontaneous imbibition it 
follows directly from Eq. ((TTJ that $x ^ H^^^ ~ t^^" and 
for the transversal width of interface fluctuations W ^ 
with [3 — x/^ where x is the roughness exponent of the 
interface (see Section llll GI for a more detailed discussion 
of X and the criticality of imbibition fronts). 

Since the velocity of the interface v — {Kpc)/{rih), the 
length scale can be rewritten in term of the capillary 
number ^x = {K/CaY^^, which highlights the importance 
of low capillary number experiments. Theoretical anal- 
ysis, to be presented in Section IIIII also shows that ^x 
is present in forced flow experiments, where it separates 
two different spatial scaling regimes. 

The fluctuations of the interface have three sources: 
variations in the permeability Sk (the fluid has to flow 
to the interface but does not do so uniformly), capillary 
forces dpc and volume (the advancing fluid has various 
pores to fill, thus Ah depends on such fluctuations). It 
will be shown in Section IIIII that disorder coming from 
saturation disorder can usually be neglected and that 
capillary and "mobility" (arising from the random per- 
meability) noise act on different length scales, charac- 
terised by ^mob ~ K^Spc/vj-jSK. On lengths / <C Cmobj 
capillary noise is dominant, while mobility disorder is 
relevant on length I ^ ^mob- Of course, the permeability 
can also be influenced by an interaction between the fluid 
and the medium and complicate this picture. 

It is also possible that for small interfacial velocities the 
physics changes. Instead of the liquid front propagation 
via jumps or advancement of the menisci, pores in contact 
may be partially wetted by film flow. From a fundamen- 
tal viewpoint this is the regime that, if described by local 
interface equations, would exhibit avalanche behaviour. 
most of the interface stays quiescent while only parts ad- 
vance. The description for such dynamics in usual kinetic 
roughening involves geometrical ideas that draw parallels 
to directed pe rcolation and its cousins in various branch- 
ing processes |216l l314j . 

III. THEORETICAL APPROACHES TO 
IMBIBITION 

The usual models, to analyse imbibition, can be di- 
vided into two main groups. There are ones built in 
close connection to experimental studies, whose goal is 
to further develop any details important for applications. 
Roughly speaking these are related to the microscopic or 
small scale properties of the imbibition systems. We shall 
also present some of these aspects in Chapter llVl when re- 
viewing experiments. Such models are particularly useful 
when the microscopic physics at the front include prewet- 
ting layers, film flow or when the saturation properties 



of the whole medium, as a function of time, need to be 
addressed. In the latter case one often needs to include 
detailed considerations about the behaviour of residual 
pockets of the gas/liquid left behind the imbibition front. 

On the other hand the models presented in this Chap- 
ter are more concerned with questions of universality in 
the morphology of the wetted region. Generally they 
include microscopic details on a more abstract ground, 
incorporated, e.g., in the noise terms included. They are 
meant to highlight the connection between the essential 
mechanisms, as ingredients into the model, and the out- 
come in macroscopic morphology. 



A. Morphological phases in imbibition 

The stability of the interface between the invading and 
the defending liquid was one of the first issues to be stud- 
ied |317| . A relation between the capillary number of 
the fluid invasion process and the size of fingering struc- 
tures in the interface has been found. In the limit of 
very thin fingers the invading fluid forms a self-similar 
fractal cluster which belo ngs to the universality class of 
invasion percolation |l54j . In fact, there is a transition 
between compact and fractal m orphology of the invad- 
ing fluid cluster |142L Il54 Il87j| . The capillary number 
Ca in Eq. Q relating viscous and capillary forces is a 
control parameter, e.g., in an experimentally obtained 
phase diagram and examples of morphologies from |l87j 
are shown in Figures |31 and ^ If the porous medium is 
easily wetted the cluster tends to be compact. The wet- 
ting properties, measured in terms of the static contact 
angle 9, are related to the fingering width of the interface. 
With increasing wetting tendency it increases and finally 
diverges at a critical contact angle 9c, below which the 
invading fluid cluster remains compact and forms a well 
defined albeit rough front |l54j . 

Compact invasion clusters with sel f-affi ne rough fronts 
have already been reported earlier |28^ . At the same 
time theoretical description and modelling as well as ex- 
perimental studies of roughening interfaces were inten- 
sively studied and have thus created a lot of interest 
in imbibition. The central objective of this Review are 
rough fronts in imbibition, and in this Chapter we present 
different model approaches to them and their theoreti- 
cal evaluation. Before that we give a phenomenological 
introduction to the elementary processes in imbibition 
fSection |lllB|l . as well as a brief recall on kinetic rough- 
ening in Section Fill CI 



B. Phenomenological approach to imbibition with 
rough fronts 

The interest on front roughening in imbibition was first 
and foremost moti vated by the assumed connection to 
the KPZ equation |147| and its underlying theoretical 
relation to experimental systems. Such spatially local 
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equations were considered after the "full" problem was 
proven to be very complicated and they hope to capture 
the essential features by focusing on the interface only. 
This description is valuable, perhaps, for (nearly) pinned 
interfaces, although also here interesting physical com- 
plications arise near the pinning transition which are the 
subject of Section lillDI 

It was however quickly realised that the detailed na- 
ture of fluid invasion could not be ignored |l2lj . The 
first theoretical efforts in front roughening with a globa l 
conservation law were done by Krug and Meakin |l6Clj |. 
They considered the roughening of a Laplacian front as 
discussed below, using a linear analy sis. This is similar 
to the problem of Saffman and Taylor |29l| , for the stable 
case of a high viscosity fluid entering a low viscosity one. 

Since the Laplacian one is a basis for later theories 
as well, it is worth considering it in some details. The 
starting point of such an analysis is based on the phe- 
nomenological law of Darcy for flow in porous media. 



VP 



(22) 



where q is the flux of liquid, rj the viscosity, and k the 
permeability of the medium, essentially dependent of the 
size of the channels through which the liquid flows (c.f. 
Equations H15(l to Ullfl in the Introduction). The incom- 
pressibility condition V • q = leads to a Laplace Equa- 
tion for the pressure 



V^P = 0. 



(23) 



The front then propagates because of mass transport by 
the flow q, whose component normal to the front deter- 
mines the front velocity 



■ 9„P 



(24) 



where 9„ represents the normal derivative to the front. 

We consider the situation described in Fig. ^1 The 
interface, of average height H, is described by a single- 
valued function y = h{x, t) with Fourier decomposition 



h{x,t) = H + ^hk{t) e 



ikx 



(25) 



At the level y = the medium is in contact with an 
outer reservoir of liquid and the pressure is supposed to 
be controlled to P(x, 0;i) = Voit). At the interface the 
boundary condition for the pressure is 



Pint = P{x, h(x, t);t) ^ Pq-p^-- (jJC 



(26) 



where pc is the coarse-grained capillary pressure, in- 
troduced in Section III Al Pq is the atmospheric pres- 
sure of the ambient gas, and /C is the curvature of the 
front. Up to linear order in the interface fluctuations 
6h{x, t) = h{x, t) — H , the pressure field is given by 

p(.,,;0=7'o-^2/-E^"^(-+-fc^^ 

K ^ V K 

fe/0 



sinh kH 



(27) 



P = Pr 








V^P = 




FIG. 15: Theoretical concepts of modelling imbibition. A 
rough front is situated at the average height H . A pressure 
field solves Laplace's equation in the bulk, boundary condi- 
tions are given at the contact hne to the hquid reservoir and 
at the rough interface. Pressure gradient is proportional to 
flux, and influx at the interface proportional to interface dis- 
placement velocity. 



where the average displacement of the interface 

^dH_^K Vq-Pq+Pc 
^ ' dt ~ ri H 



The dynamical equation for the interface fluctuation [16 
is then 



(28) 



dt 



Several important conclusions can already be drawn 
from this analysis: 

(i) The average behaviour of the interface is determined 
from the pressure at the bottom of the porous medium. A 
constant pressure gradient results in a constant average 
velocity of the interface v. On the other hand, if the 
porous medium is simply immersed into a reservoir at 
atmospheric pressure Pq, the average interface velocity 
decreases with the average height 



K Pc 

T] H' 



(30) 



Solving for dfH = v makes the Washburn-Darcy be- 
haviour apparent. 



H^{t)-H\to) = 2-pc (t-to). 
(ii) There exists an intrinsic length scale 

r]v 



(31) 



(32) 
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that separates two different modes for the relaxation of 
the interface fluctuations. If this interface moves at a 
constant velocity, this length is fixed and static, while it 
is dynamical (^x v~^l'^ ^ t^/**) for spontaneous im- 
bibition. Note that the capillary number Ca = 
relating the viscous and capillary forces, can be used to 
rewrite = (^/Ca)^/^. 

(iii) There are odd powers in momentum, more pre- 
cisely a dependence on |fc| in the dynamical Equation 
(|29|l . which reflects non- locality in space arising from 
the conservation law. However, the truly non-local 
depcndcnce cancels for wave- vectors |fc| <^ ^/H, meaning 
that the interface must be far enough from the reservoir 
of liquid to start to feel the effects of conservation. The 
surface dynamics is then local in the sense that points 
with mutual distance greater than H are essentially in- 
dependent of each other. 



(iv) This equation is in principle readily applicable to 
the case of fluid prop agation in Hele-Shaw cells filled 
with glass beads j22(Tj . but the nature of th e disorder 
remains ambiguous. Krug and Meakin 'iGd'l chose for 
their application to use non-conserved thermal noise, 
{rik{t)rik'{t')) ^ 5k-k'S{t-t'), which gives rise to asymp- 
totic logarithmic roughness, not observed experimentally 
in imbibition. Nevertheless, this equation already gives 
some hint about the influence of the conservation law. 

The simultaneous treatment of quenched disorder and 
liquid conservation was flrst attempted by Brenner and 
Ganesan |9^. They assumed that the capillary pressure 
was dependent on the spatial location x and random, 
Pc = Pc(x), as well as a constant mobility k. Their end 
result is similar to that of Eq. H29|) , although it includes 
non-linearities and the quenched nature of disorder. 



J 



dhkjt) 
dt 



= \k\ 



V H 



dk' k' [k - k'] hk-k'hk' - na dk' \k'\^ [k - k'] hk-k'hk 



(33) 



r 



The noise term r]k is obtained by assuming that the 
boundary condition Eq. (|2()(l can be modifled by adding 
a random part Pint — > -Pint + vi^i^) with so-called 
random held correlations {r}{x,h{x,t)) ri{x' , h{x' ,t')) = 
A S{x — x') S{h{x,t) — h{x',t')). These would amount to 
the fluctuating part of the capillary force at the interface. 
Note that the non-linearities also appear non-locally in 
space, contrary to the usual non-linear term introduced 
in the KPZ equation outlined in the next subsection. The 
noise term also enters the interface equation in a non- 
local way. It is however difficult to extract any results 
from this method, since Eq. H33|l is not very well suited 
even for direct numerical integration. A Flory-type anal- 
ysis (term-by-term comparison of typical magnitudes of 
the terms) yields a roughness exponent x = 3/4 but this 
result is questionable since it is well known that a simi- 
lar analysis gives wrong results in the simpler case of the 
so called quenched Edwards- Wilkinson (QEW) equation 
discussed also below. 

This anal ysis was pushed further by Paune and 
Casademunt |26lj . They considered the specific problem 
of fluid flow in a Hele-Shaw cell, where the only disorder 
is through variation in the thickness of the cell b. This 
implies that both the capillary pressure 



1 

a cos I + 



(34) 



and the permeability k cx b^/rj become random quan- 
tities. They then proceeded to derive a phenomenolog- 
ical equation that combines Eq. H33|l and earlier equa- 
tions derived by Hernandez-Machado et al. and Dube 
et al. (these are discussed in more detail later). This 



work shows that disorder in the capillary pressure and in 
the mobility act on very different length scales depend- 
ing on the velocity of the average interface, as will be 
discussed in details in Section IIII G 5|) . At the simplest 
level the main source of disorder will come from varia- 
tion of thickness b ^ b + 6b, which affects the capillary 
pressure Pc ^ Pc + 6pc and the permeability k ^ k + Sk, 



5pc 

Pc 



5r 



and 



Sk 

K 



Sr 



(35) 







K 6pc 






Pc Sk 



Based on a simple Darcy analysis, this would suggest that 
variations in the velocity due to permeability disorder 
Svk, are of the same order of magnitude as those due to 
variation in the capillary pressure Svp^ 



(36) 



which would imply that both sources of disorder are 
equally important. This may be right for the speciflc 
case of a Hele-Shaw cell with varying thickness but not 
necessarily for more general random media such as paper. 



C. Theoretical description of rough interfaces 

1. The concept of roughness 

The exact shape of a rough front or interface in a given 
dynamical process depends on the particular realization 
of randomness encountered and is unpredictable. The 
statistical properties of its shape fluctuations however 
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can be described in a controlled way. Quite often one 
observes statistical scale invariance in certain ranges of 
time and space. Naively this means that the interface 
"looks the same" when its parallel and perpendicular 
length scales are "blown up" or "shrunk" by certain fac- 
tors. In this Section, as well as in the entire Article, we 
consider interfaces with translational invariance, mean- 
ing that they "look the same" no matter which part of 
them you consider. For an illustration see Figure [TBI 

Roughness need not be a static phenomenon but 
changes its properties in time. Typically, a growing or 
moving interface will increase its roughness with time, 
as single fluctuations are accumulated. We will also see 
scale invariance in time, i.e., an interface "behaves the 
same way" when time is speeded up or slowed down with 
an according spatial rescaling. 

The relation to critical phenomena is evident. Statis- 
tical scale invariance over a certain range of length scales 
reflects itself in power law shapes of correlation functions. 
As with critical phenomena, one observes different types 
of power law behaviour of interfaces, which a re therefo re 
referred to as universality classes as well pol II6II l218j . 




FIG. 16: Different aspects of an interface, magnified by a 
factor 2 in each step from top to bottom, as indicated by the 
dashed lines between the first two panels. No inherent length 
scale becomes visible, all surface pieces "look the same" 



2. Formal description of roughness 

Let us try to formalise the intuitive picture of shrinking 
and enlarging the aspect of an interface. Generally, an 
interface S is an orientable d-dimensional object in {d + 
l)-dimensional space. A crystal surface is an example for 
d ~ 2 while a step separating two terraces on it is one for 
d — I. Since the interface is orientable, we can define an 
indicator funtion (p : M'^+^ M which takes the value 1 
on one side of the interface, —1 on the other, and right 
on top of it. In the example of a crystal one can choose 
tp = 1 inside and ip — —1 outside. 

To examine the roughness on various scales we use a 
smoothing kernel 

/i-.r(x)^^;^/(^) (37) 

for some monotonously decreasing function f{r) with the 
property dr r'^ f{r) — 1, and the d-dimcnsional 
surface volume of the {d + l)-dimcnsional unit sphere 

Bd+i ~ (rf+i)r((d+i)/2) convolution one can obtain 
via the smoothed profile function 

^i(x) ^ j d'^+'y /f+i(x-y) ^(y) (38) 

the surface smoothed to a scale L 

5L = {x|^i(x)==0}. (39) 

The profile function </5l(x) contains information about 
how far a point is away from Sl'- If 5 is a hyperplane, 
Sl remains identical to it and <pl(x) = F{5 / L) where 5 is 
the distance between x and 5, and F{r) ~ ds s*^ /(s)- 
So, it is natural to define 

6{^,Sl) = LF-\^l{^)) (40) 

as a measure for the distance of a point x from Sl ■ When 
X is taken on the original interface S Equation H4()|l gives 
an access to its roughness fluctuations on scales smaller 
than L. For an illustration see Figure [T7I 




Generally one observes power law scaling of the various 
moments of roughness, 

Wg{L,t) = (((5(x,5L)')xe5>y' - for L < m (41) 
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saturating to a constant for L larger than the so called 
correlation length £.(t). In this definition (•)xg5 takes the 
average over all points in S. The outer brackets stand 
for the ensemble average over all possible interface con- 
figurations. The measure for these averages is inferred 
from the area of the interface and is straightforward to 
construct in most cases; Xq is called roughness exponent 
of the qth moment. 

In the case of multiscaling the Xg differ from each other. 
A classical illustration of this is given by the co ncep t 
of "turbulent interfaces" introduced by Krug in |l62l |. 
Its point is best illustrated on a discrete crystal lattice: 
The height fluctuations may be related to the largest step 
heights between two neighbouring points on the discrete 
lattice. This introduces a new length scale Ah — \hx+i — 
hx\ which then makes the scaling more complicated. The 
different momenta of Wq would show a relation to the 
largest values of Aft, in L. Ah can then depending on the 
model exhibit various behaviours: E.g., its probability 
distr i buti on can be directly a function of the system size 
Il62l l303l | , or at least have a cut-off imposed by it [Til 
l235j . Many physical processes lead instead to (regular) or 
normal self-affine scaling with Xq = X — const, to which 
we restrict ourselves for most of this article, omitting the 
index q and mentioning possible cases of multi-scaling by 
remarks when necessary. 

In a roughening process the correlation length ^(t) in- 
creases with time, because physical coupling of different 
points in the interface spreads fluctuations. Again, in 
general one observes power law behaviour, which for his- 
torical reasons related to the theory of critical phenomena 
is written 



with z the dynamical exponent, 
interface fluctuations is 



(42) 



The maximal extent of 



w{t) = lim w{L,t) - ^{t)>' ~ EE t'^, 



(43) 



which connects the relations (|41|l and l|42|l and defines 
the scaling exponent /3. 

All these relations can be combined into a scaling form 



w{L,t) 



ait) atr w (^) 



(44) 



with the scaling function W{x) ^ for x < 1 and ap- 
proaches a constant for x ^ 1. In many cases the am- 
plitude a{t) is constant, which is referred to as normal 
scaling. It may however increase w ith t ime, causing one 
sort of so-called anomalous scaling '198"! . This behaviour 
occurs for forced flow imbibition with columnar disorder, 
analysed theoretically in section IIII G 61 and experimen- 
tally in Section im 

In most cases < x < 1, but it need not be restricted 
to that range. Inte rfaces wit h x > 1 sometimes are called 
"super-rough" |83 . ll98lll99j |. One-dimensional interfaces 



in disordered t wo-dimens i onal media a re examples for 
super-roughness [HI [Tsl Ell 11111113 and of particu- 
lar interest for imbibition in effectively two-dimensional 
media such as paper or thin Hele-Shaw cells [sj 



3. Height fields without overhangs 

The quantitative description of interfaces can be sim- 
plifled a good deal when overhangs can be neglected. An 
interface is then represented by a height fleld /i(x, t) mov- 
ing in time, t S M, over some d-dimensional substrate, 
whose points are x € IR''. The roughness exponent x is 
then linked to the structure factor or spatial power spec- 
trum 



5(k,i) = (|ft(k,t)|^ 



for k > l/i{t) 
for k < l/^{t) 



(45) 

where /i(k, t) denotes the spatial (d-dimensional) Fourier 
transform of ft,(x, t). The scaling behaviour of Equation 
(|45|l is illustrated schematically in Figure [T51 From this 
schematic representation the scaling form 



sik,t) ^ atf""^' s mm) 



(46) 



of the power spectrum becomes evident. As in Equation 
(gSj, we have S{k) k"^^"'' for k > 1, and S{k) = 
const in the opposite limit. 




-log^O^) -log^O^) -log^Oj) logk 



FIG. 18; Schematic plot of structure factors S{k, t) of an 
interface at times t\ < t2 < tz in double logarithmic scale. 
The correlation length S,{t), as well as the intensity of long 
wavelength fluctuations S(k 0, t) increase with time. 

Often the height difference correlation function 



G2(x,i) = {\h{x,t)~h{0,t)\' 



= J S{k, i) (1 - cos(k • x)) 



(47) 

is measured. It is self-averaging (unhke S'(k, t)), so 
"smooth curves" are easily obtained from spatial aver- 
aging, but one tends to neglect that the initial power law 
increase for |x| <C £,{t) reflects only the local roughness 
exponent xioc which due to convergen ce of the integral 
in Eq. ^ is restricted to values < 1 [US [HI. Thus, 
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to study a super-rough interface with % > 1, the power 
law decrg^se of 5'(k, t) should be analS'sed instead. 



(t ^ Lu), or the correlation function of the qth moment 
of height "jumps" over a temporal distance t 




log^Oj) log^O^) log^OJ logx 



k b 



o 




loc 



log^Oj) log^O^) log^fy logx 

FIG. 19: Schematic plot of height difference correlation func- 
tions G2{x,t) of an interface at times t\ < t2 < in double 
logarithmic scale. Part a illustrates anomalous scaling of a 
super-rough interface with x > xioc = 1, part b normal scal- 
ing with X = Xioc 

Figure 1191 contains a schematic representation of 
G2 (x, t) for a super-rough (part a) and a normally rough 
(part b) interface. In part a it becomes clear why mea- 
suring the initial power law increase does not yield the 
global roughness exponent x- 



4. Temporal interface fluctuations 

In the same way as in the spatial case, interface fluc- 
tuations can be studied in the course of time. We have 
already seen the increase of the interface width with time, 
w{t) ~ with P = x/ z in Eq. H43|) . valid for early times 
as long as ^{t) is smaller than the lateral size of the sys- 
tem containing the interface. 

Analogous to S'(k, t) and G(x, t) one can construct 
quantities reflecting the temporal fluctuation of the in- 
terface height above a substrate point x. These are con- 
veniently defined for stationary processes at late times 
in a finite system, i.e., when ^(t) = LgystQui < 00, and 
contain, e.g., the power spectrum 

5(x,c.) = (|Mx,c.)H (48) 

of the Fourier transform of the height field in time 



C,(x,t) = (|/i(x,i + s)-/i(x,s)r) 



q\i/g 



(49) 



Generally one finds an increase Cq (x, t) ~ t^^ over 
"short" time distances t. In the case of (conventional) 
scaling it is related in a simple fashion to the early time 
increase of the width in Eq. (|43|l . i.e., for all g-moments 
fiq = fi = xl z- For practical purposes Cg(x, t) is easier 
to measure than 'w{t) at early times and therefore this 
relation is often used. 

This simple picture does not hold if height fluctua- 
tions are intermittent, e.g., if the interface advances in 
avalanches with a non-trivial size distribution. In the 
experimental context this is known as Barkhausen noise, 
originally measured for the motion of boundaries between 
domains of different magnetisation in a solid with disor- 
der |2lj | . It is not surprising that imbibition into a disor- 
dered medium also exhibits multi-scaling in time, and it 
will play a role in several parts of this Review. 



5. Brief history of roughness 

Roughness of interfaces and surfaces has been studied 
in many contexts and on all scales, and the interest in 
the roughness of imbibition fronts was based on these 
previous studies to a large extent. To name a few exam- 
ples: cloud fronts, mountain rims, shore lines, rivers, for- 
est ranges, sandpiles, biofilms, cell colonies, cancer, grain 
boundaries in polycrystals, crystal surfaces, crystal ter- 
race steps, deposited thin films, flux lines in supercondu c- 
tors, fracture fronts . . . (see again e.g. pollll6lT224l314| 'l. 

The theoretical interest arose from the observation 
that many completely different roughness phenomena ex- 
hibit the same (or similar) scaling exponents, and there- 
fore be understood as belonging to the same "universal- 
ity class" . Consequently a main theoretical emphasis has 
been to flnd the essential mechanisms behind the different 
universality classes by constructing the simplest possible 
models with the same roughness properties as a given 
observed process in nature. 

To use a simple "classiflcation" there are two kinds of 
models: Continuum equations, in most cases stochastic 
partial differential equations, are more suitable to an an- 
alytical approach, either by exact solution or by renor- 
malisation and perturbative approaches. On the other 
hand, lattice or particle based models are well suited for 
simulations, often on large scales. 

Among the continuum equations we name only the 
most classic ones. First, the linear Gaussian model of 
Edwards and Wilkinson (EW equation) for a surface re- 
laxing by surface tension or downhill surface current and 
perturbed by noisy forces |89|, which can in the station- 
ary state be said to be in equilibrium. Second, intro- 
ducing a preferential direction of growth to the surface, 
Kardar, Parisi, and Zhang (KPZ) constructed a true non- 
equilibrium model. It (as they noted) is equivalent to 
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Burger's equation for an vo rtex-free, compressible, ran- 
domly driven fluid [43.ll47j . It reads 

dthix, t) = <7V^h{yi, t) + \(yh{yi, t)f -t- 7/(x, t), (50) 

where cr is a surface tension, and A measures the strength 
of the celebrated KPZ nonlinearity, proportional to the 
squared local slope of interface. In the case A = one 
obtains the EW case. ri{x, t) denotes here a simple delta- 
correlated "thermal" noise field, that roughens the inter- 
face. One has thus the correlator {r]{x,t)) r]{x',t')) = 
2T5{x-x')6{t-t'). The main ingredient of Eq. ^ is 
that the interface always has a non-zero normal com- 
ponent of the local growth velocity, with respect to the 
local interface orientation. Also, the coupling between 
the noise field rj and the non-linearity gives rise to the 
much-studied complex dynamics and morphology of the 
KPZ class. As it is in some sense the most fundamen- 
tal non-equilibrium field theory this equation has gained 
much interest and was often thought to be the correct 
large scale description for most experimental and model 
systems for surface roughening, among them for imbibi- 
tion — erroneously as we shall see in this review. 

Similarly, among the discrete models we briefly men- 
tion the family of (driven) lattice gases and equivalent 
deposition or solid-on-solid models for crystal surfaces, 
with evaporation or surface diffusion. Discrete in the 
number of deposited particles, but continuous in space 
and time are, e.g., ballistic deposition models. A peda- 
gogical overview of growth mod els c an be found in the 
lecture notes of Krug and Spohn |l6lj . Covering only the 
earlier years of studies on kinetic roughening it neverthe- 
less presents the types of models which were relevant for 
the approach of non-equilibrium statistical mechanics to 
imbibition and which are presented in this Chapter. 



D. Interfaces in disordered media and avalanches 

As already stated, there are several cases where the 
interface separates regions with common characteristics, 
distinguished by the field (j> introduced in Section Fill CI 
The example of a step separating two terraces was used, 
but interfaces will also arise between magnetic domains, 
or in paper around regions that are either burned or wet; 
in some alloys, it can be the boundary between grains 
with different orientation of the crystal structure, or even 
with completely different structure, such as in marten- 
sitic materials. 

These interfaces can move either spontaneously, under 
thermodynamic forces, or else through an applied force, 
such as a magnetic field or a pressure difference. The es- 
sential point however is the thermal influence of the envi- 
ronment is often negligible; fluctuations in the interface 
arise only from the presence of impurities or defects in the 
material structure itself. These impurities are quenched, 
i.e., static in time and the disorder felt by the interface 
is thus intrinsically dependent on its position itself. For 



example, the KPZ equation, Eq. (|50|l . is modifled to 

dthix, t) = vV^h{yi, t) + A(V/i(x, t)f + F + ■r]{x, h{x, t)), 

(51) 

where F is the force acting on the interface and the 
quenched nature of the disorder is reflected in the proper- 
ties of the noise correlator. Two broad classes are usually 
distinguished: Random Bond (RB) disorder, where the 
correlator of 77 arises essentially from the fluctuations of 
the surface tension (like in an Ising magnet with random 
coupling constants). More important is Random Field 
(RF) disorder, also since in the QEW case discussed be- 
low (Eq. (|55|l ) it can be explicitly shown that the Random 
Bond and Random Field correlators renormalise to the 
same under rescaling. 

In the RF case, the effectively felt noise depends di- 
rectly on the position of the interface and has correlations 

(?7(x,/i(x,t))?7(x',/i(x',t'))) = /\{h{x,t)-h{x' ,t'))5{x-x') 

(52) 

where A(m) is a function which is strongly peaked around 
u — Q, often taken to be a delta function 

/^{h{x, t) - h{x', t')) 6{h{x, t) - h{x',t')) (53) 

For an interface moving at some velocity v, we can sep- 
arate the average motion from the fluctuations h(x, t) = 
vt+Sh{x,t). For length scales Zth > (wi)^^-*^, the quenched 
noise correlator 

S{h{x,t)~h{x',t')) -> Siv{t^t') + 5h) -Sit-t') (54) 

V 

and the interface effectively feels an annealed thermal 
noise with effective strength A — Aq/v. 

1. Avalanches 

On length scales below Zth, the quenched nature of the 
noise is fully felt and a flnite threshold force Fc is needed 
to get the interface in motion. The velocity then shows 
critical behaviour v ^ {F ~ FcY until F ^ Fc for which 
V (x F. Equation (|51ll is then just one example of a sys- 
tem which shows a depinning transition when a control 
parameter, in this case the force F, is tuned. 

The physics in this case is best discussed in terms of 
avalanches, and it is useful to first consider the limit of 
vanishing nonlinearity A — ^ of Eq. (|51l) . the so-called 
quenched Edwards- Wilkinson equation (QEW), 

dth{x, t) = i^^^h{x, t) + F + ri(x, h{x, t)), (55) 

which is appropriate if the anisotropy (A) in the dynamics 
has a "kinematic" origin and vanishes in the limit v ~* 
0. See also the discussion in subsection IIII El about the 
origins of such behaviour and Figure QB] in Section IIVI 

The case A = has received lots of attention in the 
disordered systems community since it presents a stereo- 
type for the renormalisation group treatments of such 
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systems and associated problems. The QEW-equation is 
also an interesting candidate for a local description of liq- 
uid propagation in porous media since it contains three 
very important ingredients of the problem: surface ten- 
sion, driving force, and variations in the local (capillary) 
force, e.g., due to pore structure fluctuations. This was 
indeed noticed in the mid-eighties by Koplik and Levine 
|l57] who extended ideas of domain wa ll dy namics in 
random field Ising magnets (see also |l5l| ). And in- 
deed, we shall show below that within the length scale 
(Eq. (|32|l ). the physics of imbibition is very similar to 
the physics that comes from Eq. (|55|l . 

It is also worth keeping in mind that spontaneous imbi- 
bition exhibits something akin to a velocity-dependence: 
an interface that slows down experiences diff'erent aver- 
age velocities depending on the distance to the reservoir. 
Thus, tilts in the interface should affect the local average 
velocity, again. 

For forces just slightly larger than Fc, the interface dy- 
namics in the QEW model is characterised by avalanches. 
This implies typically that large regions of the inter- 
face stay quiescent, while only parts move in a coher- 
ent manner. This is described by the correlation length 
^11 = {F — Fc)~'^ , which is related to the fluctuations of 
the interface via a roughness exponent such that ^ Cf- 

Temporal scaling, ^ i^/^, is also assumed, and the 
length- and time scales are related via the relation 

9 = iy{z^ X), (56) 

deflning the order parameter (velocity) exponent |l92l 
|247, 249J. 

Since the disorder that the interface "experiences" has 
a history effect, it is clear that it is non-trivial to com- 
pute the correlator of the disorder in a coarse-gr ained 
QEW, as already was noticed in the beginnings |l57j |. 
Later functional renormalisation group calculations have 
resolved this issue by pushing perturbation theory up 
to second order. For uncorrelated disorder the original 
one- loop functional RG results read x = (4 — d)/3, and 
^ = 2 - (4 - d)/9, with = 1/(2 - x)- These, together 
with the ^-exponent, manifest the fact that there is only 
one temporal and one spatial scale at the critical point. 
Later work by Chauve, LeDoussal, and Wiese on much 
more tedious higher-order RG expansions [s^, Il84| have 
yielded corrected expressions for the exponents. 

For example with e = 4 — d the original C = f changes 
to C = §(1 + 0.14331e -I- . . . ). In particular in d = 1, this 
means that the RG beyond o ne loop is better abl e to ad- 
here to the numerical results [HI EH [286 ] . These 
imply that XQEW,d=i ^ 1-2 . . . 1.25, i.e., the pinned inter- 
face is super-rough. This approximate value makes for an 
interesting piece of numerology, since the exponent agrees 
with the phase held model result for the global roughness 
exponent as discussed below. 

The usual discussion of depinning in the QEW takes 
place by considering the limit F F^ and concerns the 
values of such critical exponents. One may of course ask 
questions about the behaviour for a slightly larger driving 



force, or in particular consider the avalanches right at 
the critical point. For F > Fc one has a cross-over to the 
thermal EW on large enough scales, Eq. (|54|l . and for the 
avalanches it holds that they start to overlap, eventually, 
signalling such a cross-over. It is p ossible that one can 
define here effective exponents |l57j . perhaps analogous 
to imbibition in many ways. For F < Fcina, finite system 
the velocity goes to zero, into an "absorbing state" . 

Avalanches can be observed in the vicinity of a critical 
point (for a general discuss ion about phase transitions 
with an absorbing state see [21(| , and for a n overview of 
avalanches in similar, critical systems |260j |l by imposing 
a suitable ensemble. With the caveat from above (overlap 
of avalanches) there are two set-ups that give rise to well- 
defined avalanches: 

• constant velocity ensemble, resembling forced fluid 
flow imbibition, 

• self- organised critical ensemble, induced by a com- 
bination of slowly increasing F and boundaries that 
pin the interface. 

The former is useful in this context since it exhibits trans- 
lational invariance; one may thus proceed to study in 
this ensemble the statistics of ava lanches a nd the char- 
acter of the critical point, as e.g. |l75l l248l |. SOC, self- 
organised criticality was coined originally by Bak, Tang 
and Wiesenfeld [T^ , in a context that is actually directly 
related to an interface description or model of charge- 
density waves. Recent work has expanded and clarified 
much the picture of SOC as an interface depinning tran- 
sition iSi, (^] or as an absorbing state phase transition as 
such [73 • 

The geometric picture of "avalanches" considers a 
generalised situation where the interface has "pinning 
paths" , akin to paths on the backbone of a percolation 
cluster. Assuming that there is a "punctuation event" 
if the interface is pushed at any particular spot, the 
interface is released and propagates locally till another 
connected pinning path is encountered. The geometric 
properties of avalanches set now the size distribution of 
avalanches: it is determined by the usual critical expo- 
nents and (except perhaps in the SOC ensemble) does not 
involve other physics. The important quantity is the scal- 
ing of voids, between two such paths that pin the inter- 
face, or their relation of size vs. area, which immediately 
involves the roughness exponent. The size distribution 
of avalanches is then expressed as 

P{s) 0, s-^^fis/L^). (57) 

where, if translational invariance exists (see above), D ~ 
d^ X ^nd it usually holds that Ts = 2— 2/d[^. 

In the DPD case [l35i | (see next Section and Fig. I20|l . 
it is possible to derive the avalanche size exponent in 
d = 1 -|- 1 dimensions as 

r,-l + (l/(l + x))(l-l/;^). (58) 
Similar scalings can be attempted for the SOC case. 



Note, finally, that such scaling laws as Eq. (|57|l can be 
established for the duration and support as well, not only 
for the size or area/volume. The presence of long-range 
interactions instead of the Laplacian surface tension in 
Eq. l(55|) might present some interesting twists |322j| . 



E. Cellular automata for imbibition 

Since the apparent failure of the usual kinetic roughen- 
ing framework to account for the large value of the rough- 
ening exponent was partially due to the thermal nature 
of the disorder, specific models were then developed to 
acccount for the influence of gwerac/ierf disorder. In partic- 
ular, the "Directed Percolation Depinning" (DPD) HIHl 
model completely ignores all effects of liquid conservation 
to focus only on the quenched nature of disorder. The 
algorithm that advances the front in this model is essen- 
tially similar to the sandpile (and forest fire) models Tsl 
and has three main ingredients, also illustrated in Figure 

1201 




FIG. 20: Schematic representation of the DPD model after 
0, l45| . The medium is divided into regular cells. Blocked 
cells are black, invaded (wet) cells grey, invadable (dry) cells 
white. Any dry cell adjacent to a wet one can be invaded in 
the next time step. The arrow indicates a cell which gets wet 
immediately in order to erase an overhang. 

(i) The porous medium is divided into cells that are ei- 
ther blocked or open for fluid flow. This model was first 
applied to imbibition of ink into paper, the ink pigment 
then being instrumental in blocking the fiow progression. 

(ii) The progression of the front takes place only in open 
cells, with overhangs being erased instantaneously. A 
later development of the model introduced a linear gradi- 
ent in the density of blocked cells to phenomenologically 
mimic the effects of evaporations through an increasing 
difficulty of propagation . 

(iii) The progression of the front stops when no neigh- 
bouring open cell is available, i.e., when a percolating line 
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of blocked cells is encountered. If p represents the proba- 
bility of a cell being blocked, interface pinning occurs at 
the percolating threshold Pc ~ 0.47. 

The statistical properties of the interface are related 
to the properties of the percolating path. Close to Pc the 
directed percolating path is characterised by the parallel 
and perpendicular length scales ^|| ^ \p — PcC^''^ and 
C_L \p — Pc\~'^^ , with values « 1.7 and sa 1.09. 
The saturated width of the interface then is W ^ £,± ^ 

^1^^^"" , which defines a roughness exponent x — ^^-l/^H ~ 
0.63. It was argued that the right continuum description 
of this automaton is the quenched KPZ equation, given 
above in Eq. H51|) . 

The physical justification to the main feature of the 
DPD model, the erosion of ove rhan gs has already been 
approached by experimentalists |l8nl | . It is plausible that 
propagation of the fluid in a direction parallel to the 
reservoir, if it has a chance to occur, will be favoured 
over the motion of the fluid away from it, since this does 
not result in any change in the average capillary pres- 
sure. In that sense, any overhangs will eventually be re- 
moved, although the question of the different time scales 
involved in such dynamics is certainly extremely com- 
plex. As shall be seen below, phase-field models, which 
include this possible effect do not show any kind of DPD 
roughening. Other specific predictions of the DPD or 
quenched KPZ class would be that of the DPD criti- 
cal point, a front moving with constant velocity has an 
effective roughness exponent x — 0.7 and a dynamical 
exponent z = 1. Leschhorn |l9l| has moreover shown 
numerically that the DPD model consistently gives rise 
to multi-scaling. 

A variation of the DPD idea, introduced by Snep- 
pen |306| . allows invasion always at sites of lowest re- 
sistance. An i ntere sting aspect of this model is temporal 
multi-scaling |307| , due to avalanche motion [l9d, 25^ ■ 
Avalanches are also present in the process of spontaneous 
imbibition, although the conservation law imposes a nat- 
ural cutoff on their size and distribution 81] . It is never- 
theless interesting that a similar lack of temporal scaling 
is also seen in the phase field model of imbibition, thus 
surviving the introduction of a conservation law. We re- 
fer to that in Section Fill Gl below. 

The DPD model in its simplest form may or may not 
describe well the statistical properties of a pinned inter- 
face. However, this model cannot describe the whole dy- 
namical motion of the interface, since it does not account 
for liquid conservation. It already fails by predicting a 
constant average interface velocity, contrary to Wash- 
burn's law. This also holds for dyed liquids, for which 
the DPD model was originally developed, since the con- 
servation law governing the motion of the fluid must be 
reflected on the motion of the dye particles. We conclude 
that any such local interface equation, as for instance the 
QKPZ one, is not appropriate for imbibition. Note that 
this extends to interface equations whe re the cou pling 
between the interface points is non-local |l50l l24d 1322.] . 
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F. Microscopic simulations 

The normal alternative approach to imbibition is based 
on pore network models. The porous medium is repre- 
sented by an assembly of pores and throats connecting 
pairs of pores. A simplified case is a network (in 2d) 
of horizontal and vertical capillary tubes with a stochas- 
tic distribution of cross-sections. The nodes between the 
tubes are structureless in order to avoid all the possible 
co mplic ations described in Chapter IIVI or by Lenormand 
in | l87j . With this simple stucture, the flow within each 
capillary is of Poiseuille type and the progression of the 
overall front is dictated by the global pressure difference 
and local fluid conservation at each node of the network. 

Note that there are recent lattice-Boltzmann simula- 
tions of multiphase flows that have also tackled flow in 
porous media and imbibition. The main advantage of 
the method are easy-to-implement boundary conditions. 
Other phases such as surfactant are easily added and 
questions such as Onsager transport coefficients can be 
analysed from a flrst-principle viewpoint. Actual exam- 
ples both i n 2d and 3d however highlight the technical 
challenges |20Ctl207i| associated with the method. A nat- 
ural idea is to study the front properties exactly as for 
coarse-grained models and the authours present some ex- 
amples of rough interfaces with unfortunately little anal- 
ysis that would compare to usual kinetic roughening mea- 
sures. However, it seems to us that the capacities of such 
models have to be augmented considerably before one 
can expect "realistic" results. 

A recently much practised approach is to construct 
network geometries that attempt to mimic the statis- 
tics of the porous media more faithfully via statistical 
laws for the pore sizes, and fluctuating connectivities 
(pore throats) for neighbouring pores. Other dynam- 
ical aspects (behaviour of residual gas pockets related 
to saturation, prewetting layers, Bosanquet-flow — see 
Section Hg can also be added HI ES HI IM l238l 
I279, 281, 298J. Simulations can be either quasi-static — 
the dynamics is applied at the "weakest link" , and no 
other processes are allowed — or dynamic so that the 
pressure balance is maintaine d dynamically at the single 
pore level IH US IH El [71 IHE EH ■ Such simu- 
lations are closer to describing the reality of flow through 
porous media but are naturally heavier computationally. 
In addition, quasi-static simulations have no time scale, 
which poses a genuine limitation for many comparisons 
with experiments and makes studies of kinetic roughen- 
ing impossible. Only for dynamical simulations are time- 
dependent mass flows included explicitly using Darcy's 
laws and pos sibly more complicated effects on the pore 
level IHIlS^. 

Both kinds of models have already been used exten- 
sively to study the forced and slow displacement of a 
wetting liquid by a non-wetting liquid, often in con- 
nection with oil recovery techniques. Generally there 
are difflculties for such methods to study interfacial dy- 
namics since the wetting front is a very non-compact 



object - one has to pay attention to the deflnition of 
the interface. Questions that are of relevance from 
the interface viewpoint and that have been addressed 
include pore-scale phenomena as broad pore size dis- 
tributions ,300 ] and the role of the fluid viscosity ra- 
tio of th e tw o fluids in controlling micro-fingeri ng fo r 
low Ca |335l |. studies of relative permeablities j22lj |. 
Of major interest is the use of micro-structural pa- 
rameters (from tomography or from other experimental 
sources) and the a nalys is of correlations in t he po re scale 
EllillMllMlilllMllilllMllQlll^. Typi- 
cal quantities considered are relative permeabilities and 
capillary pressures as well as saturation curves. These 
all tend to highlight the fact that in the presence of cor- 
relations the imbibition-related quantities change much 
more readily than, e.g., single-phase permeability. 

Concerning "interfacial aspects" it is of course true 
that any mass flux into the system in imbibition arises 
due to and is controlled by the interfacial area between 
the two different phases. Thus it is a natural quantity 
to compute, and consider |l08j . A re lated problem is 
the presence of initial saturations |l4lj . and some net- 
work sim ulations have b een done to address exactly this 
question flMl25l l27 ^ . 

Direct interfacial studies in network models js^, Il36| . 
Given complicated — or realistic enough — rules the 
boundary between completely non-wet (occupied by the 
non-wetting fluid) and partially saturated regions can 
be quite complex due to the effects discussed already in 
detail (film fiow, initial presence of wetting fluid etc.). 
This is of course the strength of network models. Fig- 
ure l21l illustrates the problems that ensue even in qualita- 
tive comparisons with simplistic laboratory setups. The 
generic features in the Figure show the agreement be- 
tween dynamical network simulation results and the ex- 
perimental snapshots as the fronts get more and more 
complex with decreasing Ca- Nevertheless, it is clear 
that the similarity is mostly superficial (given the many 
free choices in the details of such simulations as here). It 
would seem to be of both practical and theoretical inter- 
est to apply such network models to at least two tasks at 
hand. 

The first concerns analysing the detailed dynamics of 
the saturation in the presence of relatively well defined 
fronts — what is the actual coarse-grained driving mech- 
anism here, if the mass fiux is used as the "order pa- 
rameter"? Second, the morphology of such fronts, the 
noise affecting them and the relation to coarse-grained 
models as the phase field approach seems like an inter- 
esting avenue, in particular in the presence of pore-level 
details as snap-offs and film flow. We underline that the 
first of these suggestions has to do with the very basic 
physics of imbibition, while the second would bridge the 
gap between any analytical theory and models that con- 
tain much of the microscopic effects actually in action. 

To circumvent the difficulty of obtaining well-defined 
interfaces, conditional wetting rules, depending on the 
state of nearest and next-nearest neighbours may be in- 
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FIG. 21: Simulations of a dynamic network model vs. exper- 
imental results from a laboratory channel network |l36l Il85|| . 
The capillary numbers are Ca = 6x 10-^ Ca = 1.4 X 10-^ 
and Ca = 3 X 10"*, in the simulations, from top to bottom. 
Notice the transition to a compact front with increasing cap- 
illary number. The ratio of average pore throat conductivity 
to that of wetting layers is about 2000, controlling the "per- 
colation" of the fronts. 



troduced |31|, Il78| . These produce a more compact front 
so that a single-valued interface may be defined. This 
approach is d iscussed below, as it was used by Lam and 
Horvath |l78l |. with the goal of reprodu cing the exper- 
imental results of Horvath and Stanley |l3Cll | discussed 
in Chapter IIVI This is a priori not obvious since one of 
the main characteristics of this experiment is progression 
of the interface different from Washburn's law, v ^ H^^ 
with w 1.6 7^ 1. It was found that a very specific 
distribution of pore radii would give rise to the desired 
relation of average front height and velocity, although no 
plausible physical reason was given. In particular, it is 
not clear, that what appears to be some power law, is in 
fact transient behaviour. 

The scaling behaviour of the temporal fluctuations was 
found to be similar to the experimental results. The tem- 
poral correlation function (defined in Eq. I|49f) . here the 
moment g = 2 is considered) could be described by the 
scaling function 




FIG. 22: Schematic repr esentation of the pipe network model 
of Lam & Horvath |17^ . A regular grid of pipes is invaded 
by a fluid. Pipe diameter is random and independent for each 
pipe, determining capillary pressure at invasion and conduc- 
tivity for transport when invaded. Dark pipes are invaded, 
grey ones dry. The dashed line highlights the "interface". 



with exponents n — 0.49, 9t — 0.38, and (3 = 0.63. 

Spatially, the roughening process shows clear anoma- 
lous scaling. The stationary two-point correlation func- 
tion (cf. Eq. H47() ) had the scaling behaviour 



G2{r) = w-''5(rw(^'+'')/x) 



(60) 



where g{u) ~ for small u tending to a constant at 
larger values. The values Oi — —0.25 and x = 0.63 were 
reported. Based on the assumption that a single time 
scale controls the roughening process. Lam and Horvath 
propose the exponent identity 



(Ot + k) 



(61) 



(59) 



which, for both theory and experiment, agrees with about 
15 % error margin. 

These results raise several questions from a theoretical 
point of view. First, it is unclear why such a particular 
choice of the distribution of radii could give rise to univer- 
sal scaling behaviour. The exponent fl does not seem to 
depend in any controllable way on the disorder, so a long 
transient towards an asymptotic Washburn behaviour is 
not to be excluded. 

The value x — 0-63 for the local roughness exponent 
can be interesting since it could indicate that the imbi- 
bition process belongs to the DPD class of roughening. 
The anomalous nature of roughening and the continu- 
ously slowing down of the interface of course do not sup- 
port this view, and only a measurement of the structure 
factor could would give some further insight. 

In fact, the most questionable result is Eq. (|61|l . that 
relates the different scaling exponents. We shall show 
below how this identity conflicts with the picture of im- 
bibition that emerges from the phase field models, and 
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how the results of Lam and Horvath can be interpreted 
withm this approach. 



G. Phase-Field Models 

1. General concept 

A continuum approach different from interface equa- 
tions and conceptually related to the reasoning of Sec- 
tions nil Bl and IIII Fl was taken bvDuM et al. and 
Hernandez-Machado et al. [H H il This ap- 

proach concentrates on the bulk of the porous medium 
rather than on a hypothetical interface equation and 
presents the following advantages: 

(i) Interface equations for particular models arise natu- 
rally from the sharp interface limit of the model. Sim- 
ulations using the phase field approach thus include au- 
tomatically the relevant interface equation, with all pos- 
sible non-linearities. Equally, non-local interactions be- 
tween different parts of the interface arise naturally. 

(ii) Disorder is directly introduced at the level of the 
bulk of the porous medium and thus does not need to be 
postulated in the interface equation. 

The simplest phase field approach to imbibition starts 
from the observation that all coarse-grained descriptions 
of flow in porous media, based on Darcy's Law Eq. (|22|l . 
are essentially diffusive. All inertial effects of hydrody- 
namics are neglected and only pressure determines the 
flux of liquid q. The dynamical equation for the fluid 
fluid conc entr ation c is then defined through a continuity 
equation |295| 



dc 
dt 



V • q = 



(62) 



where the concentration is related to the saturation 
through the porosity c — sV. The concentration refers 
to both the wetting and non- wetting fluids. The two 
phases arc differentiated by the field (?i(r,t), defined over 
all space and related to the saturation concentration of 
either liquid or air as (j)^- The total concentration field 
over the porous medium is then 



c(x, t) = -c„ 



1 / 0(x^\ 



where Cw and c„ are the respective concentrations of 
the wetting and non-wetting fluids and the values = 
+4>ei—4>e) are associated to the wetting (non-wetting) 
phases. 

The field is locally conserved, which is reflected in a 
continuity equation that controls the dynamical evolution 
of the field 



9^(x,t) 
dt 



V-j = 



(64) 



where the current is taken to be proportional to the gra- 
dient of a chemical potential, j = — Af(x, 0)V/i, with a 



mobility that may be a function of both space and con- 
centration. The chemical potential has a thermodynamic 
origin, fi = SF/Scf), with the free energy functional 



F = dr 



(a(x) - ao)(j) 



and the potential 



(65) 



(66) 



For 5 = 0, this free energy represents two liquid phases, 
described by the values (/)g — r/u at chemical potential 
/i = ao. The double- well nature of the free energy, com- 
bined with the square gradient term ensures that a well 
defined interface (of width C, ^ (c/r)^/^) exists between 
the regions with cf) ~ ±(j)e- The fluids described by Eq. 
H66|l have a slight compressibility (~ t""^) but since the 
model does not consider hydrodynamical effects explic- 
itly, this represents only a minor correction to the overall 
dynamics of the interface. 

Notice that the "saturation dynamics" described by 
the phase field model is diffusive. Attempts have been 
mad e to look at th is as pect sepa rately by Mitkov et al. 
|236j (see also refs. |109| | and |237| ). This is similar to the 
approach based on Richard's equation that incorporates 
a a non-linear (and non-trivial) mobility (see Eq. 
(|2(J|l ) usually determined empirically for each particular 
case [2^. 

The model is defined in terms of the chemical poten- 
tial, so as to make close contact with theories of crit- 
ical dynamical phenomena with conserved field (Model 
B in the nom enclature of Halperin and Hohenberg |ll5j| . 
See js^ Ill2 !| for general reviews on phase fields), used 
to study binary fluids as well as the liquid-gas transi- 
tion. The correspondence between pressure p and chemi- 
cal potential /i is easily accomplished in this limit through 
p p4>e- Notice that Eq. H65|l is also a continuum limit 
of a random field Ising lattice gas with locally conserved 
magnetisation. 

Disorder enters the problem naturally though the 
transport properties of the porous medium and we can 
identify three main sources: 

(i) Capillary disorder. The term a — ao, sets the 
local value of the chemical potential and thus effectively 
represents the pressure difference due to capillary forces, 



Pc 



eOL, 



(67) 



If this quantity is taken to be random it is natural to 
assume the following correlations 



(a(x)) = a > 
(a(x)a(x')) - = (Aa)2(5(x-x'). 



(68) 



(If one assumes a Gaussian distribution for a, as one usu- 
ally does when defining a random variable by its 1st and 
2nd moment, one has to exclude negative values by hand. 
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However, the model results do not depend on the pre- 
cise shape of distribution, only on the moments quoted 
above.) 

(ii) Mobility disorder. The local mobility M(x) (ne- 
glecting any dependence on the fluid density 0), which 
controls the ease of flow of the liquid, is related to the 
permeability of a porous medium through the relation 



M 

7A2 



K 

v' 



(69) 



and is random and characterised by the correlations 
(again the precise shape of distribution does not matter) 

(A/(x)) = Ma > 
(M(x)M(x')) " A/o' = (AM)2,5(x-x'). (70) 

There is a priori no reason to expect any simple relation 
between the two fluctuating quantities and they are kept 
independent. 

Since the chemical potential is analogous to pressure, 
it also enters the model in the boundary conditions. We 
next outline the choice in two spatial dimensions, the 
only one studied so far in the literature. In the case 
of spontaneous imbibition the chemical potential obeys 
the boundary conditions /Lt(a;, y = 0) = ao (equivalent to 
atmospheric pressure) and dyii{x, y — L) — (closed at 
the upper end) js^- The boundary conditions on the 
phase field dy(j){x, y = L) — and at y = 0, the field's 
value is the solution to 



= 0. 



(71) 



Since the term ao simply corresponds to a constant shift 
in the chemical potential, we consider = from now 
on. In order to model imbibition situations at constant 
flow rate, as in |31lj |. constant flux boundary conditions 
Mdyii ~ 7o are imposed. Both situations, spontaneous 
imbibition and constant flux shall be examined in details 
below. 

When the effects of evaporation or gravity are suffi- 
ciently strong, they act in a way that stops the imbibition 
front and produce a pinned interface. Both effects can 
be taken into account in a generalised imbibition model 
of the form [sl 

= V(M(^, x)V,)-l? + ^(r, .)) 

. ^ . . . . (^2) 

The convective term is included to describe gravity. Al- 
though gravity can only be introduced properly through 
an hydrodynamical field JJO, 149J, one can argue that 
this term reproduces the correct equation of motion for 
the average position of the imbibition front. The connec- 
tion to Darcy's law is established through 



G = g -. 



(73) 



When e 7^ 0, a non-conserving term is introduced into 
the equation of motion, which, to a first approximation. 



describes evaporation of liquid (the 4> = 4>e phase), at a 
rate (jje^ and proportional to the total area covered by 
the fluid. Note, that in this case one thinks of a thin 
medium from which the fluid can evaporate on the sides, 
as in the case of a wet paper sheet. 

Taking the mobility to be M(0,x) = Mom{x.) (here 
we have neglected explicit dependence on the field value), 
the basic scales for length, energy, and time of the prob- 
lem are 



Mo^io ' 



(74) 



tm 



Equation (|72|l can then be put in a dimensionless form 
by defining 



a = — , 
Mo 

7 = G^, (75) 
which leads to a dynamical equation 

^-7^ = V(^(^, x).VMx, ^))4.(l+^(x, t)) 

^ (76) 
where the chemical potential /i = —cf) -\- (ffi — V^(/) — q;(x). 

For slow enough motion of the interface the chemi- 
cal potential is adjusted quasi-statically to the interface 
positions i which themselves are driven by the differ- 
ence between incoming and outgoing current j = — mV/i. 
In the sharp interface limit, and for a slowly moving 
front, the concentration cj) is approximately constant in 
the bulk, changing only at the interface. The chemi- 
cal potential /i thus satisfies the equivalent of a Poisson 
equation 



V • = e 



(77) 



in the bulk. At the interface, must obey the Gibbs- 
Thomson boundary condition 



(78) 



where K, is the curvature, a the effective surface tension 
of the model, the miscibility gap = — (j)- and 
Ay = V{(j)+) — V{(j)-). The quantities (j)± are the equi- 
librium values of the phase fi eld, defined by the usual 
tangent construction [381 llSlj . The interface motion is 
then determined by the normal velocity 



(A(/))ti„ = -dn^J-\± - 7- 



(79) 



These conditions are of course completely equivalent to 
the description based on capillary pressure and Darcy's 
law in the context of Eq. in Section ITlI Bl 
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2. Interface equation 

To extract the interface equation from Eq. H72I) . it is 
convenient to us e a local c oordinate system (u, s) close 
to the interface |l4Ci l332j |. Two-dimensional space is 
parametrised by x(u, s) = X(s) + u n(s), where X(s) 
is a point of the interface, n is a unit vector normal to 
the interface pointing towards the side where <^ > 0, 
and s is the arc-length coordinate along the interface 
[9l|. In terms of the phase field, this corresponds to 
^{u = 0, s) = 0. The time derivative of the field then 
becomes 
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where Vn (s) is the normal velocity of the interface at po- 
sition s. If the interface thickness, C = 1 in dimensionless 
units, is much smaller than the typical radii of curvature 
of the interface (the sharp interface limit), the Laplacian 
term of the chemical potential may be expanded such 
that 



d_ 

du 



(81) 



where /C(s) is the curvature of the interface. For a = 0, 
the one-dimensional kink solution is (f>(u^s) = 4'oiu) ~ 
tanh(M/-\/2), with corrections of order 0{pi) from the fi- 
nite compressibility of the model. Therefore, to first or- 
der, [I K. — q;(m, s) — K{s)d4>i){u) / du. Still in the sharp 
interface limit C,JC <C 1, the derivatives of the kink solu- 
tions have properties d4>Q{u) / du ~ A(j)S{u) and 



^ du 



(Po(u) 
du 



2^ (82) 



where A(j) « 2 is the miscibility gap and a is the di- 
mensionless surface tension or, for the fully dimensional 
expression. 



2 — C(Pe^J■0■ 



(83) 



Note that this quantity represents an average surface ten- 
sion associated with the interface and does not necessar- 



ily bear any relationship to the microscopic surface ten- 
sion that gives rise to capillary pressure. Nevertheless, 
assuming pc '-^ cr/ro, with tq the typical pore size and 
using Eq. yields 



ro'' 



(84) 



i.e., the dimensionless parameter a corresponds to the 
ratio between the interfacial width and the typical pore 
size. 

The capillary number C'a = Tju /a, v denoting the fluid 
velocity, represents the ratio between the viscous and 
capillary forces. Strictly speaking, this cannot be defined 
in the phase field model since the fluid is not explicitly 
taken into account, only the ratio ti/k. Nevertheless, the 
capillary number itself does not play a direct role in the 
roughening of the interface, but enters the problems only 
through the length scale ■ Assuming that the pcrme- 
alibility can be related to the average pore radius k ^ Tq , 
it is easy to show that. 



Ca 



(85) 



where v is the velocity of the interface - unlike the fluid 
velocity v. For spontaneous imbition, this is easily re- 
lated to the capillary pressure, from Eq. H84|l . to the di- 
mensionless parameter a to give Ca — v/{aa'^). 

For constant mobility m(0, x) — 1, the dynamical 
phase field equation may be inverted with the use of 
Green's function, defined by 

V^G{x, y\x', y') ^ - 5{x - x') 5{y - y'), (86) 

for the range — oo < x,x' < oo, < y,y' < oo, with 
Dirichlct boundary conditions. The half-plane Green's 
function 



l\2 



G[x, y\x ,y ) = — In — ■ — 

Att (x — x'Y -I- (y -I- y Y 
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must be used, since the presence of an reservoir at posi- 
tion y = breaks the translational symmetry in y This 
yields the equation 



/ .X' G(x|x') - i . (1 0(x', .))) . M(x, i). 

Multiplying Eq. (|55|l by \du {d(j)i^{u) / du) then effectively projects the phase field dynamics onto the interface u = 0. 
A translation u u + h(s,t) then yields 

r rh{s',t) 

ds' G{s,h{s,t)\s',h{s\t)) {Vnis')--f) + A(f>e Ids' / du' G[s,h{s,t)\s' ,u') ^ ri[x,h{x,t)) + aJC. (89) 



where ri{x, h) = Jdy 4>oiy — h{x, t)) a{x, y) ^ 2a{x, h). 



Without disorder, the advancing liquid front interface 
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remains flat, JC = 0, with position described by H{t) 
and the value of the chemical potential at the interface 
A'lint = —a, c.f. Eq. lf7H|l . It is then straightforward to 
obtain the dynamical evolution 



dH{t) 
dt 



1 



2H{t) 



7- 



eHit). 



(90) 



Without gravity or evaporation, the interface then pro- 
gresses in a typical Lucas- Washburn behaviour, H{t) ^ 
while non-zero values of 7 or e eventually stop the 
interface at an equilibrium height H.p{'-f,e) given by the 
zero of the right-hand side of Eq. (|9()|l with limiting cases 
i/p(7,e = 0) =5/(27) andi7p(7 = 0,e) = VW^- 

At early times t <^ t* = min(a7~^, e~^), the rise of the 
interface follows Washburn behaviour, H{t) — {atY^^. 
For t ^ t* , the pinning height is approached exponen- 
tially. In cases where gravity is absent, the average in- 
terface, of initial height H{t = 0) = 0, is described at all 
times by 



m 



(91) 



while for e = 0, the rise follows the transcendental equa- 
tion 



Hp J a 



(92) 



Although both gravity and evaporation pin the inter- 
face at a given height, it must be kept in mind that 



the physics behind these two effects is quite different 
|84| . Spontaneous imbibition, without external influ- 
ences, is characterised by a Laplace equation for the 
chemical potential V^/^ = in the bulk (i.e., far away 
from the interface). When solved with boundary con- 
ditions /i(y = 0) = and /i(j/ = H) = —a, this yields 
a gradient d/j, oc — in the liquid phase. The grav- 
ity term of Eq. (|72|l thus has the effect of stopping the 
interface when dfi ~ —7. 

On the other hand, the non-conserving term intro- 
duced by evaporation is such that the chemical potential 
in the bulk must be a solution of the Poisson equation 
V^/i = e and at pinning height, the gradient in the chemi- 
cal potential dfi = 0. This simply represents the fact that 
the flux of liquid from the reservoir exactly balances the 
losses due to evaporation. 

The time and length scales coming from gravity and 
evaporation can however be very different. Studies of 
capillary rise with light organic liquids (for which evapo- 
ration is presumably an extremely weak effect com pared 
to ordinary water) imbibed in filter papers |104| have 
found a pinning height of the order of 1 m, with time 
scales on the order of days. In contrast, evaporation ef- 
fects can cause the pinning height to be of the order of 
15 — 50 cm, with correspondingly much faster time scales. 

When capillary disorder is introduced, the interface 
equation is best expressed in terms of the Fourier modes 
of the interface fluctuations hk- A straightforward ex- 
pansion of Eq. (P?^ yields the linear equations 



(^hk + ie/ifc) (1 - e-^l/^lff) + \k\(H + 7) h, (l + g-^l^^l^) = i|fc| {{r^}, - aPhk) . (93) 

I 



The quenched noise enters in a somewhat non-standard 
way, by a Fourier transform of the disorder values at the 
interface, {ri{t)}k = J^e~^'^'''ri{x,h{x,t)). Immediately 
apparent is the presence of the modulus of the wave- 
vector, |fc|, arising from the conservation law. Non-linear 
terms, similar to those in Eq. (|33|) obtained by Ganesan 
and Brenner [9^ can also be obtained from Eq. (|5^ . 



3. Spontaneous imbibition without pinning 

In this section, we specialise to the case of spontaneous 
imbibition, where the porous medium is simply put into 
contact with a liquid reservoir. In this setup the interface 
follows the typical square root behaviour in time for the 
average interface height, H ~ t^^^. It is worth mention- 
ing that similar thermodynamical non-equilibrium states 
are found in other systems. An example is the invasion 



of a type II superconductor by magnetic field vortices. 
Indeed, lattice models similar in spirit to the phase field 
formalism have been developed, and it is an interesting 
question how the understanding of imbibition (dynam- 
ics, roughening, contin uum equa tions) would find appli- 
cations in that context |l38l l252j | . 

Disorder enters the model in very different ways, and 
changes in the local capillary pressure enter directly 
into the chemical potential through the Gibbs-Thomson 
boundary condition. This term is present independently 
of the velocity of the interface. On the other hand, it is 
easy to show that changes in the chemical potential dfi 
from its constant mobility value /iq due to changes in the 
mobility Sm can be estimated as 

V^iSfj.) - V((5m)V/io - V{Sm)vo, (94) 

which means that the effects of disorder in the mobility 
are essentially proportional to the interfacial velocity, a 
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quite intuitive result. 

In the case of spontaneous imbibition, the interface 
continuously slows down and the effect of mobility disor- 
der becomes smaller and smaller until finally disorder in 
the capillary forces completely controls the roughening 
process. In this Subsection, we thus consider a constant 
mobility m{(j),n) — 1, and refer to the implications of 
mobility variations in the next Subsection. 

The limit 7 = e = of Eq. (|93|l shows that there exists 
a lateral length scale separating two different modes of 
damping of the interfacial fluctuations. The changes in 
the chemical potential due to curvature (Gibbs-Thomson 
effect) and those due to slowing down of the front balance 
at a length scale 

The front is smoothed on length scales larger than 
since parts of the interface further away from the reser- 
voir have a slower velocity. Large fluctuations are con- 
stantly suppressed. 

A convenient way to study this case is to adapt the 
model to the exp erimental setup invented by Horvath 
and Stanley |130| where a steady state is reached by 
pulling a paper sheet at a constant velocity v = —vy 
towards a reservoir, or (almost) equivalently by monitor- 
ing the average interface height to remain at a constant 
value. In this case, the equation of motion reads 

dt4> = -V^ [V20 + (/} - 4>^ + a{x ~ vi)] + V • V0. (96) 

The corresponding interface remains at a fixed average 
height H = a/2v, where a freely rising interface would 
have velocity v. 

A simple spatial discretisation of Eq. H96() on a square 
lattice yields satisfactory results. The Laplace operator is 
then discretised in the standard way of finite differences 

(A<^),,, = ^ (0,,,^^.,_0,^^.). (97) 

The grid spacing Ax has to be adapted to the variation 
in the phase field and should not be larger than 1/3 of 
the interface width. 

It would be interesting to apply a finite element 
method to imbibition problems which could resolve the 
interface on a finer scale and remain coarse in the bulk 
phases. One has however to be careful to keep track of the 
quenched disorder field a{x) when rearranging the mesh. 
Also, a slightly modified version of the phase field equa- 
tions, where the disorder acts only at the interface and 
not in the bulk is better suited. This could be achieved 
by, e.g., changing the noise term to (1 — (l){x))'^ a(x) . 

For the slow diffusive motion an explicit Euler time 
integration is also satisfactory, even if the fourth spatial 
derivative requires relatively small time steps dt ~ 0.01. 
Recent developments on semi-implicit schemes, where the 
nonlinear terms are treated explicitly, but the linear ones 
implicitly, might prove useful for models of this type |88j | . 



The steady-state structure factor S{k, H) for the sta- 
tionary setup of Eq. (|^ . shown in Fig. EST a). shows 
that there is a maximum length scale for the fluctua- 
tions, independent of the lateral system size, with the 
global roughness exponent extracted from the power-law 
decay x ~ 1.25. It is slightly intriguing, that the value is 
so close to the QEW one, again (c.f. Section Fill D 1|) . 

All the data can be rescaled using ^ u^^/^ ~ 
{H/aY^'^. For the second moment height difference cor- 
relation function (defined and schematically shown in 
Fig.Hlin Section ITTrn)l one obtains the scaling function 

G2(r,F) = Aaw-x/2c,(r 1)1/2). ^gg^ 

The dependence of the interface velocity v (and equiva- 
lently on the average height H) follows from the above 
analysis. The relatively simple scaling dependence on the 
noise strength Aa is found numerically and shown in Fig. 
I23f b'). but a convincing explanation is to our knowledge 
missing ;82.. 84]. 

The scaling function g{u) is constant fo r u ^ 1 a nd 
g{u) - ff w < 1, with xioc « 1 |8ll8l ll6llT9a| . 
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FIG. 23: (a) Structure factors for systems of height H = 
50, 100, 150, 200 and lateral size L = 2H. Note that for large 
system sizes the power spectrum is not cut off by the sys- 
tem size, but by some other length scale ^x < L. (b) Cor- 
relation functions G2{r,H), scaled according to Eq. ^ for 
various v ~ H~^^^ and Aa. (c) Comparison of the correla- 
tion functions in the steady state, G2{r,H) (solid lines) for 
H = 25, 50, 100 and of the rising case G2{r, tn) (dashed lines), 
at times tn = /a. 

The same scaling picture applies in the freely ris- 
ing case, where the average interface height H{t) = 
(ai)i/2 provided that a dynamical correlation length 
Cx (t/'^y^'^ is used. There appears a complete equiv- 
alence between the interfacial fluctuations at an instan- 
taneous height H{t) and the saturated fluctuations of a 
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stationary interface. The length scale is thus concep- 
tually different from the intrinsic time dependent lateral 
correlation length commonly found in models of kinetic 
roughening (see Section IlII C|l . Here merely fixes the 
maximum range of correlated roughness and can only oc- 
cur if the "natural" dynamical exponent z < 4, so that 
the interface fluctuations can always catch up instantly 
with the available area of correlation. This also implies a 
dynamical behaviour of the total width of the freely ris- 
ing interface ^t^^'^ = t^, with /? 0.31 jll Hi- 
lt cannot be overemphasised that the exponent P found 
here is conceptually completely different from the usual 
exponent in traditional descriptions of kinetic roughening 
introduced earlier in Eq. I^Sl- 
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FIG. 24: Steady state temporal correlation functions for sys- 
tems of size Lx = 400, at heights H — 50, 100, 200. Saturation 
level and overall roughness increase with H , whilst short time 
differences and typical velocities decrease. The inset shows 
the correlation functions of moments q = 2, 4, 6 (from bot- 
tom to top) for H = 100, which decrease in their logarithmic 
slopes. 

The temporal correlations in the steady state, ob- 
tained from Eq. H96|l . can be analysed through the q^^ 
order height difference correlation functions Cq {t, H) — 
{\h{x,t + s) — /i(a;, s)!'')^/^. For small time differences t 
there appear logarithmic slopes (that have been inter- 
preted as "effective" exponents (3q) which are indepen- 
dent of H and decrease with as shown in Fig. 1241 and 
its inset. Note that the exponent (32 ~ 0.85 is funda- 
mentally different from the growth exponent associated 
with the width of a freely rising front. Standard concepts 
of multi-scaling are not easily recovered in this system, 
even if the only reason is the restriction in the system 
sizes simulated. 

It seems however certain that the reason for the de- 
crease of (3q is linked to the propagation of the front by 
avalanches, which in the QEW equation (1551) leads to 
similar statistics in temporal correlations 190]. In Fig- 
ure 1251 we show a few dozen interface configurations at 
equidistant times which apparently propagate by effec- 
tive avalanches. Zones where interfaces get blocked for 




FIG. 25: An example of an interface in the phase field model 
propagating by avalanches. Configurations at equidistant 
times are plotted together. Regions with "avalanche sweeps" 
become visible (low density of lines) , as well as pinning zones 
(high density). 



a longer time appear as thick dark lines which are sepa- 
rated by regions through which the interface sweeps in a 
fast avalanche. 

The precise nature of the decrease of jSq has not been 
studied thoroughly. An answer to the origins of such be- 
haviour may lie in the distribution functions of the fluc- 
tuating variables. In Fig. 126b ) we show the distribution 
of the normal interface velocity in the phase field model. 
It is kept in a steady state by shifting the system towards 
the reservoir at constant velocity v_^^t)ej,, as in the ex- 



periments of Horvath and Stanley [130j expressed by the 
modified phase field equation 

The velocity distribution is related to the distribu- 
tion of differences h[x,t -\- t) — h{x,t) 'm the height pro- 
file in the limit of small time steps t ^ 0. For large 
T — > oo, when both profiles are essentially uncorrelated 
the shape of the height difference distribution crosses 
over from that of Figure Fig. 126b ) to that of 126b ) . It 
seems plausible that the characteristics of this crossover 
become more apparent from the probability distributions 
themselves than from the behaviour of their moments 
{\h{x,t + T)-h{x,t)\'iY/''. The behaviour of the velocity 
distributions gives rise to analogies with other systems 
(with "local" Langevin equations) that exhibit interest- 
ing probability distributions |3^ . The deviations from 
Gaussian distributions relate to the existence of correla- 
tions, and make the question pertinent why self-averaging 
fails. 

Anomalous roughness with a global roughness expo- 
nent 1^ « 1.25 has been found in the spatial structure 
factor or power spectrum S{k,t). For the second mo- 
ments of temporal fluctuations or the exponent [32 a sim- 
ilar phenomenon is found: Figure f?71 shows the velocity 
power spectra 



S,{lo)^{\v{x,lo)\' 



(99) 



obtained from simulations of steady state distributions 
of Eq. ^ at different average heights H = 40, 80, 120, 
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FIG. 26: a) probability density function of the interface veloc- 
ity in a numerical simulation of the phase field model Eq. (|96^ . 
b) Distribution of height differences taken at two different 
times. The insets show the same data on a linear scale. 

and 160. The overbar denotes a spatial average over the 
system, the brackets and ensemble average, and v{x, uj) is 
the Fourier transform of v(x, t) = dth{x, t) which leads to 
the obvious identity S'^(a;) = uP'Suiijj)- From the power 
spectra one can read a global temporal exponent /3|'°''^': 
A possible power law behaviour goes like 

(100) 

It cannot be overemphasised that this is the global coun- 
terpart to the local exponent defined via Cq^2it, H) in 
Figure I24L and not the exponent related to the early 
time increase in interface roughness when starting from 
an initially flat configuration. 

The crossover observed in the velocity distributions, 
Fig. 1261 will of course prevent a "pure" power law in 
Sy{uj), but we can nevertheless observe three regions giv- 
ing rise to "effective" exponents, (i) For low frequencies, 
long times, the fluctuations are stationary, Sv{uj) ^ uj 
and therefore f32 = 0. (ii) Fluctuations at short times 
and high frequencies, on the time scale o f ava lanche du- 
ration, are "turbulent" in the sense of |l62j |: We find 
Syiu) - and thus pf"''"' = 5/4. It certainly 

would be of interest to examine these in greater detail 
and pin down the connection to the velocity distribu- 
tions shown in Fig. 1261 (iii) Last we can identify an in- 
termediate regime with an apparent behaviour between 
Sv{uj) ^ uj~^^^ and cj^^/^, which would correspond to 
^global _ g^g^ would however suggest not 
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FIG. 27: Power spectra for the interface velocities in the 
steady state of Eq. 19611 . At low frequencies Sv{i^) ~ t^i, at 
high frequencies S4uj) = ^^1-2/3!'°'"" with f3f°^''^ « x ~ 5/4 
indicating a dynamical exponent z ~ 1. The behaviour in 
the intermediate frequency range is not clear. Data were ob- 
tained in systems with size L = 500 at average interface heights 
H = 40, 80, 120, and 160 respectively 



to put too much emphasis on the power law nature of 
this regime, which probably is caused by an interplay of 
several different mechanisms during the crossover from 
avalanche propagation to the stationary fluctuations at 
long time scales. 

The different regimes appearing in Sy{uj) exhibit differ- 
ent scaling behaviour with respect to the average inter- 
face velocity |v|. The stationary regime with 
ends at ~ |vp2 . A possible interpretation is to 
take the typical duration of avalanches Tava — A/iava/^ava 
from their vertical extent A/iava ^ their "sweep- 

ing velocity" WavaS?m|v|~''2 The velocity Uava also 

appears in the large frequency tail where Sy{uj) ~ uj~^^'^ 
with an amplitude oc v"^^^. The lower bound of this tail, 
LUi , at the crossover to the intermediate regime seems to 
be independent of the average velocity |v|. Also these 
findings certainly need to be examined more thoroughly. 

Interestingly we find the exponent ^|'°'"'' = 5/4 in the 
short time regime to be the same as the global roughness 
exponent x- Assuming that these short time fluctuations 
follow a standard scaling picture we can combine them 
to z = x/beta ~ 1 which means linear or ballistic prop- 
agation of fluctuations along the interface. At first sight 
this is plausible: We are dealing with the advancement 
of the front in the course of single avalanches, between 
encountering regions of stronger pinning. Here the front 
can propagate freely with some typical velocity, and it 
is easy to imagine that fluctuations will be carried along 
the interface with a constant velocity during that pro- 
cess. Again, this would be an interesting topic for further 
studies. We have not pursued an analogous direction, by 
computing the power spectrum of the interface velocity 
(recall Fig. l26|l . In the context of microscopic models - re- 
lated to e.g. Barkhausen noise - with avalanches that are 
easier to define than in imbibition, the S{uJv^int) can be 
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related to the exponents describing the avalanches which 
in turn relate t o tho se that describe the underlying crit- 
ical point itself |l67j | . 

It is worth remembering here both with respect to the 
correlation functions Cq and the /3-exponents, that the 
deviations from ordinary kinetic roughening stem from 
the memory or history effects in imbibition. In local mod- 
els the central question in this respect concerns whether 
the local fluctuation 5h{x, t) = h{x,t) — h{t) defines effec- 
tively a Markovian process. Imbibition is different, since, 
among others, clearly the dynamics measured by the two- 
point functions Cq depends on both the fluctuating field 
Sh{x, t) and also on h, the distance to the reservoir. Sec- 
ondly, it is clear though only on a qualitative level that 
the existence of implies memory effects that are valid 
also beyond the particular details of spontaneous imbibi- 
tion. 

In any case, the model results shown in Figure [^re- 
produce seve ral ex perimental features found by Horvath 
and Stanley |l3Cl| |: First, the late time saturation level 
of Cq increases with H indicating larger overall rough- 
ness, and second, Cq{t, H) at fixed small t decreases with 
H, indicating a faster intrinsic avalanche velocity for a 
smaller H . 

However, the picture developed so far contradicts 
strongly with the exponent identity, Eq. proposed 
by Lam and Horvath |85| . This identity is derived un- 
der the assumption that the dynamics are controlled by 
a single time scale which can be obtained by two means. 
For a moving interface, a width w is reached after a 
time ti ^ On the other hand, the form 

of the time correlation function, Eq. (|59|) suggests a time 
scale t2 ~ ~ y^(et+K)/«:/3^ ^^leie 9t = 0.38 and 

/3 = 0.63 are found. Assuming that ii cx ^2, the exponent 
identity is obtained. This is trivially true in usual kinetic 
roughening but wrong for spontaneous imbibition since 
ti and t2 describe two distinct physical time scales. It 
takes a time ti for the correlation length to have value 
£,v{ti). This in turns controls the width as w ^ On 
the other hand, t2 is the relaxation time of the fluctu- 
ations when the interface is kept at a fixed height. In 
this case, the correlation length is a predetermined 
constant and the saturation within this zone is obtained 
when the spatial extent of the correlations equals ^« . We 
emphasise that the growth of the width in spontaneous 
imbibition is controlled by the increase of , the zone 
available for correlated fluctuations, and not from the 
intrinsic dynamical fluctuations of the interface. 

In fact, the correlation function, Eq. H60() . clearly de- 
fines a length scale ^„ ~ , where 7 = {9i + K)/a ~ 0.4, 
which then implies G2{r) ~ ^ig{r/£,v)- Thus, ^„ ~ u^^ "^ 
is analogous to the ^x ~ v~^^'^ discussed for the ideal 
case of spontaneous imbibition with Darcy- Washburn be- 
haviour. Equation H6()(l also defines a global roughness 
exponent, x = ko:/{k + Oi) = 1.25, again very similar to 
the results obtained from the phase field model. 

The scenario proposed by Lam and Horvath could only 
occur if the fluctuations were slower than the increase 



of Cx, which implies t2 ^ ii- Recall the point that 
for z < 4, in the case of the phase field model, this is 
never true. Indeed, the numerical data indicate ^2 ^ w"^'^ 
and ti ~ nf^-"^^ which shows this assumption to be false 
asymptotically. 

4^. Spontaneous imbibition with pinning 

In presence of the pinning effects of gravity or evapora- 
tion, numerical simulations of the phase field model show 
that the behaviour of the average interface height is well 
described by the transcendental equation, Eq. (|91|l and 
Eq. (|92|l . The actual pinning heights deviate only slightly 
from the predictions due to the presence of disorder |85j . 

It is convenient to analyse the effects of gravity and 
evaporation separately. In the absence of evaporation, 
the linearised equation of the fluctuations Eq. 193|l . im- 
mediately shows the existence of a length scale (t) re- 
stricting the range of the spatial fluctuations and evolving 
in time as 

where H{t) is given by the solution of Eq. (|92|l (minor 
corrections of order 0{£,g/H) are also expected). The 
length scale (,g{t) is thus analogous to the length scale 
^x (t) in spontaneous imbibition without external influ- 
ences, but it does not increase in time according to a 
simple power law any more. At pinning, this length scale 
becomes 

4(^p) = ^(^)'^'=ex(i/p). (102) 

As in spontaneous imbibition without gravity, dynami- 
cal scaling relations can be established by assuming a sin- 
gle correlation length £,g{t) ^ 7^1/2 [aH/Hp)^^'^ , where, 
from Eq. (|92|l . the quantity H/Hp is a function of j'^t 
only. The temporal behaviour of the width of the inter- 
face W{t) can then be scaled as 

W{t) = -f-^/^w{j^t) (103) 

where w{x) ~ x^^'^ for x 1 and tends to a constant 
for large arguments. Numerically, the global roughness 
exponent x ~ 1.25 has the same value as for spontaneous 
imbibition without external influences js^ . 

The scaling picture at pinning is confirmed from the 
structure factor. Figure [2HI shows the structure factor 
for systems with different pinning heights (arising from 
different values of the constant 7). The structure factor 
decays as k~^'^ at large values of wave- vector, consistent 
with the value x = 1-25 and the curves can be collapsed 
on the common scaling form by setting £^„{Hp) ~ 7~i/2 ^ 

Hi/': 

Sik,Hp) = (eg(i/p))i+2^s(fceg(i?p)). (104) 
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The scaling function s(fc^) is the same as in the freely 
rising case and thus relates to the scaling function g in 
Eq. 

As for pure spontaneous imbibition the concept of a 
dynamical exponent z describing the propagation of fluc- 
tuations along the interface does not apply to imbibition 
with gravity. The fluctuations always catch up to the 
available zone of correlated roughness, correlation lengths 
cannot exceed ^g. 

Dynamical roughening of the interface in the presence 
of evaporation (and absence of gravity) is the most com- 
plex case. Although the average interface dynamics is 
well described by the mean-field equation, the extra term 
in e(l — e^^l'''!^) means that a single correlation length 
^ = ^{t,H{t)) cannot be unambiguously identified. 

The situation becomes clear only in the pinned inter- 
face limit, where a simple correlation length can be de- 
fined as 
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aC' = 2e(l-e- 



)• 



(105) 



In the limit ^ Hp, ~ i'^/^Y^^ while in the opposite 
limit ^ Hp, ~ (a"^ /{ae)y/^. The first case corre- 
sponds to weak evaporation, defined by e ^ c? ja^ , while 
the second case is that of strong evaporation, e ^ 6? ju"^ . 
The value of the roughness exponents is established from 
the decay of the structure factor, with again the result 
X = 1.25 113.^ 

The prediction for the pinned correlation length can 
be checked by considering the scaling behaviour of the 
pinned structure factors 5(fc, i7p) as shown in Fig. 1^ 
The data with the highest evaporation rate can all be 
collapsed assuming a correlation length {oLt)"^!'^, 
while for lower evaporation ~ e^^l'^ . The crossover 
occurs roughly for parameters a = 0.2 and e = 10~^. 
The complete set of data can be reduced to a common 
scaling form by solving Eq. (|105|l numerically (using a 
value of cr = 2-\/2/3). As a result the structure factor at 
pinning has a scaling form similar to Eq. Ij 104(1 , although 
care must be taken in identifying the correlation length. 

To summarise, the main result of the analysis pre- 
sented here is that the fluctuations of the interface are 
intimately coupled to the average position, through the 
length scale • Both and H should be studied simul- 
taneously. Experimentally, the existence and dynamical 
behaviour of this length scale can be inferred from the 
saturation of the spatial and temporal correlation func- 
tions and should be easier to observe than the precise 
values of the scaling exponents on scales r < ^x ■ 



5. Constant flow 

Imbibition with a constant flow rate can easily be per- 
formed, e.g., in Hele-Shaw cells with disorder and can 
be modelled within the phase field approach by a sim- 
ple modification of the boundary conditions. So far 
we have had for the average velocity of the interface 
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FIG. 28: A log-log plot of structure factors of the pinned in- 
terface in presence of gravity. The inset shows the value of 
S{k, Hp) for pinning heights Hp = 20, 40, 60 and 80. The 
main figure shows the collapse of the data under the assump- 
tion of a correlation length (,g{Hp) ~ H^^ . The dashed line 
indicates a roughness exponent of x = 1.25. All quantities 
are taken in the dimensionless units of Eq. H72|l . 
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FIG. 29: A log-log plot of structure factors of the pinned 
interface in the presence of evaporation. The main figure 
shows the value of S{k, Hp) for pinning heights (from bottom 
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to top) : Hp ~ 20 {a = 0.2, e = 10"''); Hp 
0.2, e = 5 X 10-*); Hp ~ 45 (a = 0.1, e = lO""*); Hp ~ 63 
(a = 0.2, e = 10"*); Hp ~ 77 (a = 0.3, e = 10"*); Hp ~ 109 
{a = 0.3, e = 5 X 10~^). The structure factors for e ~ 10"" all 
collapse on the same curve, independent of the value of a. A 
complete collapse of the data can be accomplished by solving 
Eq. 11051 numerically, as shown in the inset. The dashed 
line indicates a roughness exponent of x = 1-25. Again, all 
quantities in the dimensionless units of Eq. 17211 . 
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2v = |(V/^|y=o)| (the factor 2 coming from the misci- 
bihty gap). If now the disorder is in the mobiUty, its 
effect on the interface dynamics is proportional to the 
flow influx, so a priori it cannot be neglected. This adds 
a term to the interface equation, of the form 



j dx' G(x|x') (V • m(x)VAi) . 



(106) 



The usual projection technique becomes quite involved 
in this case. However, for a linear analysis, it is sufficient 
to relate the gradient of the chemical potential to the 
average interface velocity V/i ^ v, which yields 

hk ^-{v + afc2) \k\hk + i |fc| {rj}k + V {m}k (107) 

where {m}k ~ 2 J_^e^^''^m{x,h{x,t)) represents the 
noise induced by the mobility. In contrast to the cap- 
illary disorder, represented by |A:|{ry}fc, this noise is non- 
conserved, but with a strength proportional to the mean 
velocity of the interface. 

Th is equation, also derived in l26l| and partly in 
|l22j | shows that mobility and capillary disorder influ- 
ence the fluctuations on different length scales. Small 
length scales are dominated by capillary disorder while 
the effect of mobility disorder becomes dominant on large 
scales. This is evident and natural, since the latter re- 
flects fluctuations in the local fluid flow in the bulk, which 
obviously have an effect beyond local pore randomness. 
A crucial point however is that the crossover length-scale 
fmob between the two regimes is velocity dependent, in 
terms of dimensional quantities 



smob ^ r- 
VT] OK 



(108) 



In the simple model of Paune and Casademunt, where 
both permeability and capillary pressure are related to 
an effective pore radius ro, fmob ^ ro/Ca |26lj . 

The length scale fmob can be seen in Figure 1201 where 
the roughness spectra of three model setups are plotted 
on top of each other: (i) capillary disorder only, which 
gives the main contribution to roughness on short scales 
< Cmob, (ii) mobility disorder only, which roughens the 
opposite regime, and (iii) including both, which gives an 
envelope of the other two curves. The comparison be- 
tween the structure factors corresponding to different av- 
erage interface velocity is shown on Fig|^ An effective 
exponent /3 ~ 0.4 can be defined, but it is strongly influ- 
enced by cross-over between the capillary and mobility 
regimes of disorder. 

The behaviour of the interfacial fluctuations is a lot 
clearer at low capillary number, where only capillary 
disorder is relevant . The intrinsic length scale = 
{a/vy/^ is still present but it is now static, and does not 
play an explicit role in the dynamics of the fluctuations. 
Figures and show the numerical results obtained 
from a simulation of the phase fleld model with constant 
flow, rn(x) = 1 and without pinning (7 = e = 0). The 
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FIG. 30: Comparison of roughness spectra (structure factors) 
for capillary and mobility disorder. The dashed line represents 
capillary disorder only, the dot-dashed line mobility disorder 
only and the full line to both sources of disorder present. 
The simulation is done at i; = 0.005, taken at an average 
height H — 50. Beyond a length scale ^mob (cf. Eq. IIO8I1 ) 
mobility disorder, with x ~ 1 gives the dominant contribution 
to roughness, on smaller scales disorder in capillarity, with 
X = 1.25 is more important. 



structure factor at saturation shows that the short wave- 
length exponent x ~ 1-25 with a much smoother rough- 
ness at large length scales. The dynamical behaviour 
of the width shows flrst a steep increase, with exponent 
(3 ~ 0.42 which is indicative of a dynamical exponent 
z — xl P = 3. At later times, the width flattens out, to 
an effective exponent (only deflned on the small window 
available) (i « 0.1, more representative of logarithmic 
roughening than true power-law behaviour, a trend con- 
firmed by the form of the structure function. 

The dynamics is thus controlled by two correlation 
lengths, ^1 ^vt, which dominates the dynamics for wave- 
vectors fc <C and ^3 ~ (crt)^/'^, which controls the 
scale k ^ The time development of the width is 

determined by the ^3 term as long as ^3 ^ while the 
slower ^1 is only effective at times when ^1 ^ ^x ^ 



w{t) 
w{t) 



tx/3 



t > ty 



(109) 



The crossover time between the two dynamical regimes 
is tx ~ cr^/^/u'^/^ and the width after saturation of the 
rough zone w{tx) ^ ■^x • 



6. Columnar Disorder 

Recent experiments have also d ealt with the case of 
imbibition with columnar disorder [mill, defined by 

(a(x)a(x')) - = {AafS{x-x'). (110) 
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FIG. 31: Comparison structure factors for both capillary and 
mobility disorder at different velocities. The simulations are 
done for three different velocities, u = 0.005, 0.01 and 0.05 with 
the structure factors taken at the same average height H — 50. 
A clear separation between the two roughening regimes oc- 
curs only if ^mob ~ i/y oo. Beyond a length scale ^mob 
(cf. Eq. 110811 ) mobility disorder, with x ~ 1 gives the domi- 
nant contribution to roughness, on smaUer scales disorder in 
capillarity, with x ~ 1-25 is more important. 
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FIG. 32: Growth of the width for interfaces moving at con- 
stant velocity v = 0.005, 0.01 and 0.05. The inset shows an 
initial growth exponent j3 ~ 0.42 followed by a weak loga- 
rithmic regime. The main figure shows the rescaling of the 
data using Eq. H109|l and a crossover time tx ~ v~^'^. The 
roughness exponent x = 1.25. 



An example is shown in Figure 15^ fSection llV A l|l . The 
main characteristic of this type of noise is that the di- 
mension perpendicular to the reservoir does not come 
into play, it corresponds to imbibition on "stripes" of 
different properties. Thus, at the level of the interface 
equations, the position of the interface does not enter 



FIG. 33: Structure factor at saturation for interfaces mov- 
ing at constant velocity v — 0.005, 0.01 and 0.05 taken in the 
logarithmic regime for the growth of the width (nearly sat- 
urated regime). The inset of the Figure shows the structure 
factors, the faster velocity have the smaller range of corre- 
lated roughness. The main Figure shows the rescaled data 
using the correlation length ~ v^^^'^ and the roughness 
exponent x = 1.25. 



the quenched noise. The linear interface equation is thus 



dt 



= ^{v + ae) \k\hk + - \k\Tjk (111) 



where now {rjkrjk') = {Aa)^Sk,-k' ■ With this simple na- 
ture of the disorder, the roughening process can be stud- 
ied analytically at the linear level. At saturation, the 
structure factor 



Sk 



1 



(w -I- crfc2)2 



(112) 



indicates a short- wavelength roughness exponent x = 3/2 
and a flat interface at large distances. However, due to 
the development of fingers in the interface, separated by 
large slopes, the linear approximation of Eq. (|lll|l breaks 
down and experiments (Figure l¥T)l . as well as simulations 
of the phase field model (Figure ISljl . show instead intrin- 
sic anomalous roughening 198]. In this case, there is a 
time-dependent increasing amplitude of the structure fac- 
tor, characterised by an exponent 9a, due to steepening 
of local slopes on the interface. 
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-2(x+e„)-i forfc»r' 
for k < 



(113) 



and ^ — t^/^. Indeed, the e xperimen tal distribution of 
slopes shows a power law tail |31lll312| . similar to what is 
observed for roughening in fractal structures [Of. At this 
point it becomes clear that the linear interface equation 
fails to appropriately describe this case. 
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FIG. 34; Structure factor obtained from the phase field model 
with columnar capillary disorde r an d constant flow rate. In- 
trinsic anomalous roughening |19^ is visible in the inset. 
The scaling plot (main graph) indicates a roughness exponent 
X ~ 1.25 and dynamical exponent 2 « 2. The scaling scaling 
plot indicates an anomalous exponent da ~ 0.75, similar to 
the value found experimentally |312ll . The simulations were 
done with a imposed velocity v = 6.25 x 10"* and capillary 
disorder only. 



IV. EXPERIMENTS ON IMBIBITION 

In the following we summarise the current experimen- 
tal situation in imbibition with regard to the theoretical 
ideas presented earlier. The main goal of this section is 
to give some experimental substance to theoretical ideas 
presented earlier. Beyond this we want to understand 
where further interesting avenues exist, and to outline 
which of the recent advances in models have been con- 
firmed. From the existing works one can see some clear 
directions for further investigations. 

The introduction presented already some experimen- 
tal complications associated with both spontaneous and 
forced imbibition. Since our approach is predominantly 
based on an interfacial approach, where the advancing 
front of liquid has an effective surface tension, it is im- 
portant to discuss in which cases the simplifications of 
the theoretical models are expected to be valid. 

To recapitulate the viewpoint of a statistical physi- 
cist interested in the front roughening, the main issue is 
whether an interface description actually makes sense at 
all. In our case this is considered to be true if the large- 
scale geometry of the boundary between the "dry" and 
"wetted" parts of the system can be described with the 
aid of self-afRne concepts. Exceptions can arise if the sat- 
uration of the imbibing fluid is a very smooth function, 
in the coarse-grained sense. Then the transport dynam- 
ics are determined on the level of individual pores, and 
coarse-graining to an interface is not possible. When an 
interfacial description is possible, pore-scale phenomena 
can be embedded in the noise contained in the interfacial 



equation of motion. A second question to be asked is 
whether the Darcy- Washburn behaviour is valid (recall 
that it is based on a idealised picture of fluid flow), even 
in the presence of a well-defined imbibition front. Failure 
of simple Darcy-flow concepts can stem from fluid-matrix 
interactions, from hydrodynamical reasons (inertial and 
transient effects, as Bosanquet-flow), and from contact- 
line physics, li ke the change of the contact angle over 
time |2Ml3n^ . 

In this section, when appropriate, we also point out 
analogies between kinetic roughenin g in imbib ition and 
other experimental systems, (see |2y, lllfil l224j for early 
reviews of the subject). We first discuss the specific 
problem of paper, or fibre networks, and then continue 
with results on interface dynamics for constant flow rate 
(forced fluid flow) imbibition, followed with the sponta- 
neous case. 

Before embarking on this task, it is worth keeping in 
mind some basic facts of roughening in nature. The the- 
oretical models are always a "coarse-grained" description 
of the experimental system considered, that is the micro- 
structure enters only via an effective description, most 
often encoded in the noise terms. Several length scales 
are also usually present and many cross-over phenomena 
can occur. In real systems, as well as in simulations, 
one encounters a typical structural length-scale Ic (e.g., 
the bead size in a Hele-Shaw cell) so that real scaling 
behaviour is established only for I ^ Ic- Cross-over be- 
haviours are also due to the finite size L of system it- 
self which enters the scaling function of any measure of 
roughness or correlations. Finally the boundary condi- 
tions can induce boundary effects by suppressing or en- 
hancing local dynamics. It is customary to neglect these 
by considering only a part of the system (to use a window 
in space), but there is of course no guarantee that this 
procedure works a priori. In particular, it is often as- 
sumed that one needs to subtract the 'tilt' (linear trend) 
from the height profile h{x, t) to avoid long-wavelength 
effects (see e.g. 297] for a discussion). For all these three 
reasons it is clear that one should take any exponent val- 
ues measured for a fixed l/lc with a grain of salt. Indeed, 
the effect of boundaries on the behaviour of systems ex- 
hibiting criticality is a very complex topic by itself, with 
many facets, e.g., for d irect ed percolation (183j and the 
KPZ universality class |l63l |. 

A second important point is that true scaling be- 
haviour is observed onl y if good averaging, in the "ther- 
mal" sense is obtained |245j . The interface must see sev- 
eral different realisations of disorder before any scaling 
behaviour can be observed. This is intuitively clear even 
in a simulation, where one run just considers a particu- 
lar noise history and representative scaling is established 
only within a statistical ensemble. There are two ways 
around this problem. The simplest, brute force approach, 
is simply to do a large enough number of experimental 
runs with different noise configurations. This is how- 
ever difficult with imbibition since, analogous to fracture 
experiments, it is often impossible to repeat the experi- 
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ment very many times. The other option is to consider 
long enough experimental runs, so that th e interface fluc- 
tuations reach a steady state |l3fl l245j . This is again 
very difficult in imbibition, due to the slow time scales 
involved. Also, in the theoretically interesting case of 
spontaneous imbibition the front slows down constantly 
and there is no steady state. In cases where thermal 
noise becomes irrelevant (slow, avalanche- or activation- 
like dynamics) the situation is in practical terms even 
more precarious than in usual kinetic roughening when 
it comes to obtain reliable data. 



7. Difficulties inherent to paper imbibition: a case study 

In many statistical physics experiments on imbibition 
the medium of choice has been either a Hele-Shaw cell 
filled with (transparent) beads up to a certain volume 
fraction or, for obvious simple practical reasons, paper. 
The first approach is perhaps less subject to accidental 
experimental complications. To both assist with the in- 
terpretation of actual scaling results and to highlight the 
difficulties that could arise in various media we now dis- 
cuss the specific case of paper and its impact on studies 
of imbibition. Sheets have been used to study striped 
media |35lj | and radial p ropa gation 225J, and two fronts 
can be made to collide |246l |. in addition to the papers 
cited below in later sections. 

Experiments have demonstrated that clear, Washburn- 
like behaviour can only arise if there is not much inter- 
action (swelling, prewettin g) b etween the fibre networks 
and the penetrating liquid |l04j . Ordinary paper is made 
from wood fibres and, in many cases, chemical additives 
and filler material s like talc and clay, arranged in a dis- 
ordered structure '253'|. There is a wide distribution of 
pore sizes (this can be simulated by compu ter m odels 
|254j . or studied via, e.g., X-ray tomography \29^ ) with 
a high effective tortuosity. The actual mass distribution 
per un it ar ea varies (due to a phenomenon called floccu- 
lation |253| ) and the resulting structure has on the top 
of this power-law correlations up to a cut-off of a few 
average fibre lengths |269j . 

The surface of paper sheets is extremely uneven and 
the concept of surface pore is itself ill-defined. This gives 
rise to problems in defining static and dynamical con- 
tact angles as is typical of inhomogeneous substrates and 
implies that there will be large fluctuations in the cap- 
illary pressures. In commercial papers, there is also an 
anisotropy between the "machine-direction" (MD, par- 
allel to the web) and "cross-direction" (CD, perpendic- 
ular to the web), which may also influence the liquid 
penetration in the structure. For industrial papers there 
are often residues of chemicals, which may induce drastic 
changes in the effective viscosity or surface tension of the 
invading liquid. 

In paper, the whole fibre structure may be modified 
by swelling of the fibres in contact with the liquid [s^ . 
Cellulose fibres show a great affinity to water and can 



absorb large quantities in millisecond time-scales, giv- 
ing rise to concomitant changes in fibre volume and pore 
structure. By nuclear magnetic resonance (NMR) it is 
possible to observe the simultaneous intake of water into 
pores between fibres but also inside the fibres themselves. 
However, this is not the case for many organic fluids and 
oils. Any intuitive picture of fibres (or the network) as 
capillary tubes is thus often false: the pore structure is 
high ly non-tri vial and in some cases time-dependent (see 
also IllllMI). If swelling occurs, the volume to be filled 
with liquid increases and the flow resistance of the pore 
network changes, which leads possibly to non- Washburn 
behaviour. 

There can also be an exchange of liquid between the in- 
side of a flbre and the "surface" pores, important in par- 
ticular since the imbibition takes place "along" the sheet. 
This complicates enormously the fluid flow since there are 
no well defined "structures" (either the pores, or the fi- 
bres) responsible for the capillary forces ^12. J5^ .263j . 
The fluid may very well undergo imbibition as such only 
on the rough paper surface. This is a mechanism that has 
bee n studied itself by various authors in other contexts 
mini but not for such rough surfaces. AU this 

makes it questionable as to whether it is possi ble to eas - 
ily find a representative network description |l86l ISOOj l . 
Mathematically, the problems amount to defining the es- 
sential effective properties of a pore for the imbibition 
process. Classical applications of imbibition to study the 
fiow properties of porous media hope to answer this un- 
equivocally, but correlations in the pore-scale properties 
may however render the starting point useless, since both 
the typical scale of pores and the effective correlations 
need to be ma thematica lly specified in any model of a 
"typical" pore p^llsf)^ . 

A. Forced flow imbibition 

Simultaneously to the development of the theoretical 
picture of "local" kinetic roughening several groups at- 
tempted to measure the prop erties of a fluid-gas interface 
in Hele-Shaw cells [ol Il2ll l288j. These were typically 
done in planar cells of thickness between one and a few 
bead diameters. The general outcome is that no really 
convincing link could be established with such "local" de- 
scriptions. Quite recently this has been taken up by the 
group of Ortfn, Hernandez-Machado, and co-workers, us- 
ing specially prepar ed substrates with controlled disorder 
|l22ll31ll I312L l313l | . We outline below the current exper- 
imental situation. 



1. Interface roughening 

The early works ^ EH measured typically the 
temporal exponent /3 and the roughness exponent x by 
using digitised pictures of interfaces and two-point cor- 
relation functions (c.f. Eqs. 14711 and H49|) '). This method 



has many natural limitations: the typical resolution 
available for spatial structure was of the order of 250 
points of h{x, t) and the two-point correlation function 
cannot measure roughness exponents larger than 1. Rep- 
resentative values for the exponents would be /? = 0.65 
and X ~ 0.80 . . . 0.9. These are substantially higher than 
standard K PZ e xponents /3kpz = 1/3 and xkpz — 1/2. 

He et al. |l2ll | observed the roughness in a Hele-Shaw 
cell with capillary numbers 10^^ < Ca < 10"'^. They 
determined a large roughness exponent x 0-8 — 0.9 
(as illustrated in Fig. 133 from 121]) and a decrease of 
the global width with increasing capillary number with 
a possible saturat ion o f the width at very low capillary 
number (Fig. OEI jl2l|). The roughness exponent was 
obtained from the two-point correlation function and is 
probably indicative of a super- rough (x > 1) interface. 
They also made the important point that the dynamics 
is non-local and introduced a simple scaling argument 
to relate the total width of the interface to the capillary 
number Ca 

w ^ C-" (114) 

where q = 2^/(2% -I- 1) in 2D. This argument is ob- 
tained by balancing the local (pore level) capillary forces 
to the global viscous pressure drop driving the interface 
and agrees qualitatively with the experimental results at 
large capillary numbers, although it cannot fit with the 
plateau at low capillary numbers. 

However, this argument does not take into account 
any large scale curvature of the interface, from which 
the length scale ^ Ca in Eq. 1)21(1 arises. The ex- 
perimental results can then also be explained with the 
theoretical picture developed in the preceding section. 
Since the experimental capillary number is small, we as- 
sume that only capillary disorder is important and that 
the mobility has a negligible effect. At very small cap- 
illary numbers, the length of the cell L <C and so 
the width is independent of Ca. As the capillary number 
increases, < L and, neglecting the small logarithmic 

correction at large scales, w ^ ^ Ca^^^ ■ A roughness 
exponent 1.0 < x < 1-25 yields 0.5 < 0.625 which is also 
in qualitative agreement with the experimental results. 

At large capillary numbers, Horvath et al. noted that 
there seemed to be a cross-over in x to values ~ 0.5 
|l28j at large length scales (a clear example is depicted 
in Fig. 1371 from the same reference) which would be con- 
sistent with a regime where either mobility disorder or 
the "thermalisation" of the quenched disorder becomes 
relevant. The measured dynamical exponent f3 — 0.65 is 
also indicative of a smaller value of z. 

Lately the forced flow case has been studied intensively 
both theoretically (see the preceding section) and experi- 
mentally b y the group of Ortin , Hernandez-Machado and 
co-workers [HI Iml El 113. The central idea, differ- 
ent from earlier approaches, was to use Hele-Shaw cells 
with a pattern of copper tracks or islands placed on a sub- 
strate of fibre glass. An illustration of the experiment is 
shown in Fig. 1381 together with the way of introducing 



36 




CAPILLARY NUMBER C 



FIG. 35: Interface width vs. the capillary number Ca in a 
Hele- Shaw cell with a constant flow rate. (Experiment by He 
et al. |l2ll ]'l. The lines show scaling predictions (see text). 
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FIG. 36: Two examples of th e loc al width (called W here) for 
the experiment of He et al. |l2lj . For L small the effective 
roughness exponent is about 0.8, whereas for larger scales the 
scaling changes. 



randomness into the samples. Here such disorder is given 
simply by the fluctuations in the cell gap. A definite ad- 
vantage of this setup is that the disorder is under better 
qualitative control. In particular it allowed for the first 
time the study of "columnar disorder", where the gap 
value is quenched and does not vary in the direction of 
interface propagation, for the first time (see Eq. (|110|) '). 
On the other hand this very interesting work also high- 
lights the standard experimental difficulties in this field: 
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FIG. 37: The spatial two-point correlation f unct ion in a forced 
fluid flow experiment by Horvath et al. Il28|l . Notice the 
presence of two scales in the data, with a larger exponent 
at small scales. 



the need for many samples (disorder averaging) and the 
problems in attaining reasonable scahng windows for the 
possibly critical quantities are difficult to circumvent. 
The in practice achievable length-scales are obvious in 
Fig. 1391 and one should note in particular the case with 
columnar disorder. 

Hele-Shaw cells with a pattern of square islands 
(square disorder) is similar to using cells with randomly 
packed beads. Intuitively one may perhaps expect a bet- 
ter qualitative picture of the porosity (gap) fluctuations 
than in such earlier setups. For capillary numbers in the 
range 1.33 < Ca < 17.0, the cro ss-ov er phenomenon al- 
ready observed by Horvath et al. |l28l |. together with ex- 
ponent (3 ~ 0.5 was also present. The structure factor of 
Fig. 1401 clearlv shows a short-range exponent XSR '^-^ 1.0. 
At large scales, it appears that the roughness exponent 
changes continuously as a function of velocity, but this 
behaviour should be contrasted to the theoretical re- 
sults of Section IIIII From the phase-field simulations 
(see Figs. QUI and 1^ . it is known that roughening due 
to capillary forces gives rise to a large roughness expo- 
nent but only up to a length scale • At larger length 
scales, mobility disorder dominates, but there can be a 
very long crossover between the two regimes, particularly 
at large capillary numbers (i.e., for a large interface ve- 
locity). The experimental data of Fig. ^Ula-re then more 
indicative of a well-defined long-range roughness expo- 
nent xlr ~ 1.0 and of a cross-over behaviour, than a 
true velocity-dependent roughness exponent. 

For columnar disorder, analysed theoretically in Sec- 
tion IIII G 51 the papers brought up the idea of a phase 
diagram, which depends on the strength of the capillary 
forces (gap in the cell in the experiments being the con- 
trol parameter) and the forced flow velocity. If the cap- 
illary for ces are of major importance, anomalous scaling 
is found |l98l l312| . This results from another length- 
scale, the "lag" between neighbouring columns due to the 
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FIG. 38: Sketch of the experimental setup of 3111. a) and b) 
present views of the experimental apparatus, c) a SEM image 
of the copper islands on the fibre-glass plate, and d) and e) 
present views of the disorder pattern for "point-disorder" and 
"columnar disorder" , such that copper islands are denoted by 
the black regions. 



difference in local mobility. The usual signs of anoma- 
lous behaviour could be extracted: correlation functions 
for the spatial fluctuation would exhibit different scal- 
ing exponents Xq depending on the order (or moment) 
of the correlations measured, the local and global ex- 
ponents for the roughness would differ, and finally the 
global width would in addition exhibit "super-rough" be- 
haviour, X > 1 198]. The scaling of the structure factor 
of the interfaces also pointed out the possibility of sev- 
eral X values, as shown in Fig. l41l for the experiments: for 
short scales the roughness exponent would exhibit such 
anomalous, super-rough scaling {x ~ 1-2) while on more 
macroscopic scales the same cross-over behaviour as for 
square disorder would be observed. Again, this may be 
a crossover effect between capillary and mobility disor- 
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FIG. 39: Interface profiles for examples of disorder, point- 
like in a), columnar in b), with the original images included. 
Notice in particular the scale of the vertical "discontinuities" 
in b). From |31l|l . 
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FIG. 40: Results from Ref. |31l| for the power-spectra of the 
interfaces at different times approaching a steady state under 
point disorder as in Fig. I39f a'l . 



der. For the other exponents /3 w 0.5 and z « 2.0 were 
obtained (sec Fig. 142(1 . Again, higher order correlation 
functions show different effective exponents, a possible 
indicator for multi-scaling. Fig.l43lgives evidence for this, 
for interfaces with columnar disorder and at saturation. 
A distinction between the various phases for multi-scaling 
is illustrated in Fig. ^| Such signs are similar to what is 
seen in the phase-field description of interface dynamics 
(though the actual exponent values may not agree), and 
indicates that the presence of non-local dynamics invali- 
dates simple scaling. 
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FIG. 4f : Results from |31lj for the power-spectra of the inter- 
faces at four different imposed velocities and columnar disor- 
der. The curves are shifted for clarity. The vertical line gives 
the scale corresponding to the typical width of the columnar 
tracks, and the other lines are fits to determine roughness 
exponents. Notice how a large velocity seems to imply in par- 
ticular a larger short-range roughness exponent (recall that 
the structure factor decays as A;"'^"'"^''-'). 



2. Avalanches and depinning 

All the early works in particular in Hele-Shaw cells 
pointed out more or less directly the difficulty in observ- 
ing the slow, avalanche- like motion for small Ca, due to 
the lack of self-averaging. Further experiments then con- 
centrated on the study of avalanches and noise properties. 
Horyath et al. studied the noise distribution P{rj{x,t)) 
|129| . For local interface equations^ this can be obtained 
directly from the local increments of the interface height 
{v(x,t) ~ ri{x,t)), as can easily be seen from, e.g., the 
slope-averaged KPZ equation. This quantity, turned out 
to have a power-law tail which is in qualitative agreement 
with the fact that no standard exponents were observed 
either. Fig. 1451 shows the power-law-like shape of P{rj)). 

In later papers by Dougherty and co-workers, the 
avalanche-like motion of the interface was studied in more 
detail 0, |^ by looking at the average interface velocity 
as a function of the local average slope |V/i| (inside a 
measurement window). It would be first assumed that, 
for a fixed average interface velocity, the function {v) 
is parabolic in the interface gradient. This is a conse- 
quence of the non-linear, KPZ-like terms in the inter- 
face equation which imply (u(A/i)) cx A(A/i)^ in the 
ensemble-averaged sense. As the velocity of the inter- 
face approaches zero, depinning models predict that the 
quantity (i;(A/i)) can either vanish, which indicates the 
isotropic depinning universality class (QEW-like), or not, 
meaning that the depinning is anisotropic (QKPZ-like) if 
one considers local interface equations. The parabolic de- 
pendence of the velocity on the slope was indeed obtained 
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FIG. 42: Results from ,311. | on the roughness exponents, 
where the type of disorder is varied (circles: point-like, trian- 
gles: columnar). The top one shows variation with imposed 
velocity, and the bottom one with the gap between the plates 
between which the advancing liquid is confined. Solid and 
open markers correspond to short and long range exponents, 
respe ctive ly. For x > 1 (or o in the notation of the figure, 
&om |31l| [') denoting "super-roughness" the region is shaded 
in grey. 
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FIG. 43: q-th order two-point correlation functions of the in- 
terfaces in the experiment of Soriano et al. show at saturation 
some type of effective multiscaling as their logarithmic slope 
(taken as Xq) is clearly q-de pendent. Here the disorder used 
was of columnar type |312(l . The long-range scaling is not 
reliable due to poor statistics. 
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FIG. 44: An experimental phase diagr am for anomalous scal- 
ing in the setup of Soriano et al. |312)l . The symbols indicate 
the points explored. The numbers next to the symbols tell 
the difference between the "global" and "local" roughness ex- 
ponents, X Bt.nd xioc. The effect is strongest for slow interface 
progression and small values of the gap — which renders the 
disorder strong and would perhaps indicate a simultaneous 
decrease of the effective interface tension. 



experimentally, and the results further indicated that the 
depinning process was isotropic 6], consistent with the 
value X — 1.25 obtained from the phase field simulations 
(see Fig. 

Meanwhile in another paper, using the same setup, it 
was established that the avalanche statistics are not quite 
akin to normal depinning models: the area distribution 
was close to an exponential one and not a power-law, 
which also was manifest in the fact that the avalanches 
spread mostly laterally |8l| . One constructs from images 
of interfaces (Fig. I47|l typical avalanches by using time 
steps (Fig. I48|l to recover the areas, locally, swept by the 
interface. Clearly, the data in Fig. 0^ are more in line 
with an exponential distribution. 

This is in clear distinction from critical behaviour 
found in self-organised systems (e.g., SOC sandpiles or 
QEW). Criticality in these systems is measured in terms 
of the order parameter (average velocity) going to zero. 
An interface is then expected to exhibit avalanche-like 
motion with a wide distribution of scales both perpen- 
dicular and parallel to the orientation of the interface. 
Another landmark would be a power-law distribution 
of avalanche sizes that also is a signature of the lack 
of an intrinsic scale. These experimental results imply 
at least indirectly that the non-locality of the interface 
dynamics, which provides a definite length scale for the 
avalanche dynamics, was of importance, as one would ex- 
pect in particular in the limit u — > 0. They also highlight 
the difference to other scenarios (theoretical and exper- 
imental) where the details of microscopic physics upon 
coa rse-g raining yield a "local" (KPZ) interface equation 

iliiliii. 



In this respect we also would like to point out the possi- 
bility of measuring directly (in model media) quantities 
that relate to the non-local properties as pressure flue- 
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FIG. 45: The effective noise distribution P(rf), or equivalently 
the distribution of the local interface velocities in a 2d forced 
fluid flow imbibition (open circles). Clearly, there is a wide 
range of local rates of advancement. For comparison the data 
include the result from the RSOS model (stars), tha t bel ongs 
to the KPZ universality class (see text, figure from |128| '). 



tuations and saturation |331j . These could be perhaps 
combined with various analysis of the interface dynamics 
in future experiments. 
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B. Spontaneous imbibition 

So far, both the theoretical and experimental devel- 
opments in imbibition highlight the specific features of 
fluid invasion together with quenched roughness, and it 
is important to outline what in particular one should look 
for in the case of an interface that spontaneously slows 
down. Some of the specific ideas are: 

• a length-scale, as , separating regimes of effec- 
tively local and non-local interface dynamics, 

• multi-scaling, or anomalous scaling, absence of sim- 
ple scale invariance, 

• the interface velocity and its variations as the in- 
terface approaches an equilibrium height, 

• scaling functions for fluctuations, including the 
properties of pinned interfaces. 

Next we overview both the older experiments and the 
more recent ones, while trying to relate them to the back- 
ground established. 



FIG. 46: a) Dependence of the local velocity u{v, s) on the 
local slope s. The average velocity is given in bead diame- 
ters/second. In b) the same data is represented after rescal- 
ing by the velocities The parabolic shape of u{v, s) would 
imply the presence of a KPZ-like nonlinearity in the effective 
dynamics. 




FIG. 47: Typical interface in a forced fluid flow experiment 
by Dougherty and Carle IsH . The velocity was 13.7 /im/s, 
which should be compared with the scale in the figure. 




1. Pinned interfaces 

A set of statistical physics imbibition experiments were 
done by Buldyrev et al. |3|4l,|43 and Family et al. [95l |. 

The first one was performed with a dye solution in a ver- 
tical capillary rise setup. The rising front moved from the 



FIG. 48: The contrasted image (from 0) shows where the 
interface has moved during a 4 s interval. As can be deduced, 
the motion takes place rather uniformly in small patches, 
without the presence of a wide variety of sizes of "avalanches" . 
This might be taken to indicate the presence of a correlation 
length, due to the non-locality of the dynamics. 
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FIG. 49: The distribution of avalanche sizes in the experiment 
of ref. 81] . The tail resembles more an exponential than a 
power- law. 



reservoir, until eventually, the dye front stopped, due to 
gravity and/or evaporation (no dynamical measurements 
were done) , and the roughness of the pinned interface was 
measured. The main experimental finding was a rough- 
ness exponent x = 0.63, consistent with the DPD/KPZ 
model. It should be noted, however, that the length scale 
of the scaling behaviour was extremely small: For a of to- 
tal lateral extent L = 40 cm, C{1) ~ only for length 
scales I smaller than Zmax ~ 1 cm, with a cross-over to 
a constant or to logarithmic behaviour at larger /. This 
scaling region is only a few times larger than the average 
fibre length in paper, so the result should be taken with 
some care. A similar experiment was done in a three 
dimensional sponge- like material [4^ . The stopped in- 
terface yielded a roughness exponent x^^"^^ ~ 0-5; again 
consistent with the 2d DPD model. 

Amaral et al. V| studied the role of evaporation in- 
duced pinning, by controlling the pinning height via the 
evaporation rate (presumably by modifying the humidity 
during the experiment). The result is that the width of 
the pinned interface is related to the pinning height hp 
(and thus to the evaporation strength) through a novel 
exponent 7, as Wsat ~ h^- The experimental value was 
found to be 7 = 0.49, which could also be related to 
a modified version of the DPD model as shown in 
Fig. [Sn| together with the scaling ansatz used. This re- 
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FIG. 50: The scaling of the interface width for a pinned in- 
terface as a function of the evaporation rate, related to the 
equilibrium height of the interface. Both the "raw data" and 
the scaling used by Amaral et al. are depicted 0j. 



suit can also be related to the length scale £,x{h): if 
the role of evaporation is to stop the interface, with lit- 
tle or no infiuence on the statistical fluctuations, then 
^W~e^(M'-/*P^'and7 = x/2. 



Kumar and Juma |166| also performed an experiment 
in presence of varying evaporation conditions, with the 
claim that the roughness exponent depended strongly on 
the evaporation rate, i.e., x = xi^)- 

The last results in this respect are those by the Mexi- 
can group 0,^3- Here, the pinned interfaces in paper- 
ink wetting experiments were analysed. The main idea 
was to compare the front roughness (after usual digitisa- 
tion procedures from a digital picture) between samples 
of different size. The pinned interface width vs. sample 
width was stated [ig'I to follow an exponent of the order 
of 0.75 or 0.63 . . . 0.64 depending on the paper grade 
(type of sample). It is worth noting that the Washburn- 
behaviour was absent, i.e., the measured interface heights 
vs. time increased much slower. In such conditions it is 
no surprise that the interfaces possess lots of overhangs. 
Such a result is, in any case, qualitatively in accordance 
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with the fact that Washburn-scahng is not observed - this 
must imply, that the pinning height can, as again was ob- 
served, vary with the sample size as well. However the 
later interpretation ,17] that the actual exponent value 
would relate to the correlations in porosity and thus on 
the local "force" acting on the interface must be ques- 
tioned. Unlike in many rocks etc. the pore space in pa- 
per is constrained by the finite thickness (of the order 
of 0.1 mm) and a typical pore thus has a size of a few 
tens of /X m^, while the pore space is not expected to 
have fractal properties on any sizeable scale. There is 
no evide nce t hat pores are correlated in paper. (Note 
that ref. |269l | discusses ureal mass fluctuations and the 
correlations there in, not pores.) 

These results highlight the fact that the slow dynamics 
close to pinning are very difficult to understand. Indeed, 
Delker, Pengra and Wong, and later Lago and Araujo, 
have demonstrated in Hele-Shaw cells that the approach 
to pinning can follow a power-law, similarly to depinning 
transitions at the critical point in general Il7fit Il77| . 
This would follow after a Washburn-like temporal regime 
has been observed. With pinning allowed or included, 
Washburnian dynamics would imply an exponential law, 
also in contradiction with such experiments. 

Recall that at the critical point of an absorbing state 
phase transition, or a depinning transition, the veloc- 
ity often follows a power-law decay in time. The actual 
exponents of the experiments are perhaps not very sur- 
prisingly close to any standard depinning models. The 
attained values found by Delker and Wong for the decay 
exponent of the interface velocity were in excess of unity, 
and dependent on many factors. Indeed, similar power- 
law behaviour has been reproduced for a contact line in- 
side a capillary tube. The Hele-Shaw one may of course 
also simply demonstrate how local evaporation or slow 
surface flows lead to creep-like motion of the interface 
|47| . Other experiments on spontaneous imbibition in 
such systems ha ve addres sed saturation and cross-overs 
to fingering flow |l20l l282l | . These results do not seem to 
be in line with any of the standard local models, which 
does not seem to be too surprising, either. 



2. Moving interfaces 

In analogy with the pinned interfaces, the results on 
moving ones are also quite varied. Family et al. per- 
formed an experiment in a horizontal capillary setup 
with water, with the position of the air-water interface 
recorded both temporally and spatially 95,]. The main 
results were an average interface progression h ^ , 
with S = 0.7 (faster than the Washburn one) and a self- 
afhne interface described by a Family- Vicsek scaling re- 
lation, and characterised by the exponents /3 = 0.38 and 
X = 0.76. Concerning h{t) ~ it may be so that the 
details of the setup were at play: a reservoir was placed 
underneath the paper used, to combat evaporation. This 
may have caused condensation, or prewetting, thus in- 
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FIG. 51: The temporal two-point correlation function from 
the stationary height experiment of Horvath and Stanley, 
|l30l . The apparent roughening exponent is about 0.56 re- 
gardless of the speed {V = 4.03 x 10"^, 1.32 x 10"^ and 
1.40 X 10"'^ cm/sec). Saturation takes place the sooner the 
larger the velocity. 

creasing the velocity of the front. It is of course no news 
that water penetration in paper can be non-Washburn- 
like. The spatial scaling regime was rather small, for 
distances below Imax ~ 2 cm for a 40 cm wide sheet of 
paper. 

The temporal scaling of the interface was studied in 
detail by Horvath and Stanley |130( . A paper sheet was 
moved so as to keep the interface always at a fixed dis- 
tance h above a reservoir. A power law behaviour for 
the time correlation function C2{t) was established, with 
^ , (3 = 0.56. The experimental result is shown in 
Fig. 1511 The velocity v at which the paper must be moved 
towards the reservoir varied as t; '--^ h~^, — 1.6. This 
implies that the interface propagated slower than t-^^^, 
h{t) ^ ti/(^2+i)^ from V = dh/dt. No spatial scaling re- 
sults were reported. The value of (5 was independent of 
h, in contrast to the saturation time of C2{t). The cor- 
relation function was brought into a scaling picture, 

C2{t) - w-®^/(to(®'+®^)/^'), (115) 

where f{y) is a scaling function such that /(y) ~ 
for y <C 1 and f{y) ~ const for y ^ 1. The values of 
the independent exponents were Ql = 0.48 and 8t = 
0.37. One can again try to rewrite this in terms of by 
relating the lateral length to the velocity at ft,, so that 

which yields C2(t oo) ^ h^/"^ ^ v^^l"^ . As long as 
^x < the correlation function Ciify is independent of 
the total width of the system. Horvath and Stanley used 
only paper of one, fixed, size. 

A similar velocity result was obtained by Kwon et al. 
|170( . A paper towel was set on an inclined glass plate. 
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to follow the capillary rise of a dye solution. This gave 
h{t) ~ Two scaling regimes were present in the sat- 

urated width; on small length scales (< 2 cm), x = 0.67 
while on larger length scales (up to 20 cm) x « 0.2. 
Within the simple Family- Vicsek picture they obtained 
f3 — 0.24 on the short length-scale regime. It is conceiv- 
able again that the cross-over is simply due to non-local 
effects. 

The only experiment so far on horizontal front rough- 
ening was performed by Zik et al. (35Q| |. They ob- 
tained rough interfaces only with what was called "highly 
anisotropic" paper, with x — 0-4- For isotropic paper, 
the roughness was at best logarithmic. It is remarkable, 
that the scaling for the anisotropic paper was observed 
through a large range of length scales, not only up to a 
few fibre lengths. 

There are two recent, intriguing experiments on spon- 
taneous imbibition that have shown promise with respect 
to agreement with theory. First, Soriano et al. also used 
the Hele-Shaw setup to explore the effects of w = by 
preparing the interface by finite-velocity transients. This 
means that the interface was first forced, as in a fixed 
flow rate experiment, until a certain height, and then 
let to relax to v = |312| . It was established, with 
columnar disorder, that the local progression inside the 
copper tracks (with preferential wetting properties) fol- 
lowed Washburn-like dynamics. This setup and, possi- 
bly, experiments on regular assemblies of tracks would 
certainly be useful to investigate the cross-over between 
local and non-local dynamics (again, ) and the role of 
annealed (thermal) disorder by repeating the test in the 
same fixed geometry, as illustrated in Fig. 1521 where the 
effect becomes visible in the local width w{l,t) and the 
saturation power-spectrum. 

In a recent paper Geromichalos et al. explicitly demon- 
strate the existence of a scale-length similar to fx by 
studying rising liquid fronts between two rough glass 
plates fiOlj : we reproduce examples in Fig. First, 
the Washburn-scaling of the velocity was established, al- 
though with standard prefactor-problems: the dynamic 
contact angle is unknown quantitatively. Measurements 
of the spatial two-point correlation function indicated the 
presence of two kinds of scaling, with XSR ~ 0.8 and 
Xlr ~ 0.6 (see Figure IS^ . Maybe coincidentally, the 
values and the existence of the cross-over are similar to 
what has been established i n a number of forced-flow ex- 
periments discussed above |l2ll Il28l l312j | . 

The cross-over scaling should follow, as a function of 
the cell gap d the law 
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which was in line with the results for dodecanol and wa- 
ter as the rising liquid. The exponents differ from most 
of the models except for perhaps that of Ganesan and 
Brenner |99| . and also one should note that xlr is nu- 
merically close to the quenched KPZ value of 0.63. Unfor- 
tunately the published results do not consider the other 
exponents, as e.g. (3 (unpublished data may imply, that 
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FIG. 52: Interfaces with a regular spacing of the tracks (or 
disorder columns) 6 mm wide, a) snapshots at different times, 
b) local width w{l,t) for two values of I, and the power spec- 
trum at saturation. Th e vertical line indicates wavelength of 
the pattern. From |312(| . 



P ~ 0.5 . . .0.6 ■1031) nor the structure factor or higher 
order correlation functions that could be used to study 
multi-scaling. The authors pointed out that fronts some- 
times had clearly visible dynamical structures or "kinks" 
one of which is shown in Figure 1551 This would be cer- 
tainly worth further study as would be the repeatability 
of the measurements and correlations between runs using 
the same glass plates several times. 



C. Imbibition, wicking, and liquid transport in 
living organisms 

Transport of water and solubles, and water entering 
parts of an organism are essential processes of its life. We 
finish this Section by giving a (very) short overview on 
recent studies of imbibition and water transport in plants. 
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FIG. 53: A rising spontaneous imbibition front in an experi- 
ment with two rough, pa rallel glass plates. The snapshots are 
separated by 10 seconds llC 



focussing on those which by their methods or subject 
come close to the physics aspects, as presented in this 
Review. 

In a plant water flows through capillaries in the stem, 
roots, branches, and leaves, in vessels forming a fractal 
structure. Perhaps not too surprising from a Statistical 
Physics point of view, scaling behaviour can be observed 
in this structure for various quantities, among them the 
vessel diameters in relation to the biomass, the ratios of 
biomasses between roots and leaves, between stem and 
leaves, etc. Some simple arguments indicate that the 
plant structure, as it is found, o ptim ises the supply of 
cells with respect to resource use |9^ . 

One has to be careful to apply concepts of physics to 
biological systems, e.g., Darcy's law connecting pressure 
gradient to flux, because living organisms may, and gen- 
erally do actively influence all kinds of processes — so 
the exchange of water and solubles between a cell and its 
surroundings. Interestingly, the term imbibition is used 
for the soaking of water into seeds and grains, a mainly 
passive process and therefore in line with the models and 
experiments presented above. Imbibition to seeds starts 
germination through different mechanisms. Anybody is 
familiar with the swelling of a wet grain, which breaks the 
outer hull and thus enables or eases the seedling's growth. 
But water entering the grain also starts physiological and 
developmental processes nec essary for the early develop- 
ment of the embryo [1 IH [HI ^E^. It is influenced 
among other paramet ers crucially by temperature and 
water supply |35l l344j . Surface properties influence wet- 
tability and imbibition Besides these agricultural 
issues knowledge about grain imbibition h as great irnpor - 
tance in food processing and food storage Il39l l290j | . 

Nuclear magnetic resonance allows for a direct mea- 
surement of water flow inside the grain. Generally, the 
hydration of a seed is a multistage process. For example 
in soy beans the water concentration keeps increasing on 
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FIG. 54: In a) the spatial correlation functions of the exper- 
iments of Geromichalos et al. .lOl] are demonstrated, with 
varying separations of the two glass plates. The presence of 
two scaling regimes is apparent, as is a cross-over length-scale 
that decreases with the separation. In b) the length-scale is 
shown, both water and undecanol, as a function of the gap 
width. The square-root scaling of Eq. (11171 is roughly repro- 
duced by the data. 



the time scale of days with water entering different com- 
partments one after anot her and reaching the embryonic 
axis relatively late |264j . Similar result s are found for 
cerea ls in de taile d studies, so for barley |llOl f223]. oats 
[isH . maize |22^ . 

It is of considerable interest for agricultural applica- 
tions to understand the interaction of a seed and its sur- 
rounding soil during imbibition. It would allow for op- 
timal preparation of the soil structure (grain size distri- 
bution) and its humidity and give an estimate measure 
for successful germination. Brucklcr has built a model 
based on a Darcy ansatz for imbibition of maize grains 
and was able to determine essential parameters such as, 
e.g., grain wall permeability by direct measurements of 
sing le grains which were kept in well defined conditions 
|40| . Subsequent tests in seed beds yielded good corre- 
spondence of observations and model predictions under 
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FIG. 55: A detailed look of some fronts (with undecanol), 
showing the presence of coherent featur es or deterministically 
moving "kinks" in the front dynamics ([Toi|). 

variation of soil moisture and soil grain size, but difficul- 
ties in estimating the role of temperature variation |34| . 
This is in fact not too surprising, as temperature has sub- 
stantial influence on physiological properties. As stated 
above, active processes in living organisms involve many 
details, often too complicated to be captured by simple 
thermodynamics or generally physical reasoning. 

V. CONCLUSIONS 
A. Overview 

In this Article we have reviewed various aspects of front 
roughening in imbibition. To us, with a statistical physics 
viewpoint, it is striking that imbibition appears in many 
different fields, from the most abstract statistical physics 
of rough interfaces to detailed engineering studies in, e.g., 
paper fabrication or oil recovery. Man y yea rs have passed 
since the original work of Washburn [239'| , but perhaps 
it can be stated that the general understanding of imbi- 
bition has advanced relatively little. 

This work has aimed at collecting the knowledge ac- 
cumulated so far, grouping it into fields, and stating the 
most urgent and interesting open questions and advances 
from the theorist's viewpoint. By far the most work has 
been done in numerous detailed technical applications, 
as the references to engineering literature on the subject 
that we have gathered should bear witness to. 

Imbibition has attracted a good deal of interest among 
statistical physicists, in great part due to the presence of 
experiments on interface roughening that (still) remain 
unexplained. As usual in this context, it is of particular 
interest to look for the simplest possible model that cap- 
tures the essential ingredients of the problem. This how- 
ever implies the presence of the usual dichotomy between 
theory and experiment: The simplified approaches of sta- 



tistical mechanics predict scaling behaviour but have dif- 
ficulties to make quantitative predictions and there is of- 
ten no satisfactory connection between experiment and 
theory. This is true in part for imbibition and it remains 
to be seen how quickly the distance between theory and 
experiment can be overcome. This also holds for the 
comparisons between numerical simulation models that 
contain microscopic details - as outlined in the Introduc- 
tion - and more coarse-grained ones like the phase field 
approach. We next summarise the main theoretical is- 
sues that have been resolved so far and then point to a 
number of topics and issues that remain open and deserve 
further attention. 

Theorists, in evaluating experiments and model sim- 
ulations, have quite uncritically been looking for power 
laws in surface roughness without looking at the global 
picture. In our opinion, this has in many cases obfus- 
cated the view for the essential features of imbibition. 
For example, it is of primary importance to observe and 
quantify the lateral length-scale (cf. Eqns. and 
(|95|l '). coming from the interplay of surface tension and 
pressure gradient in the liquid bulk. We note that this 
comes out of theoretical considerations js^, Ei H| , and 
that there is recent experimental pro of for its existence 
in the case of spontaneous imbibition |10 ij , while forced 
fluid flow in Hele-Shaw cells yielded indications even ear- 
lier ITp . 

Imbibition problems are of a broad general theoret- 
ical interest since they present a field in which usual 
"non-local" interface models (here referring to equations 
of type dth{x,t) = f_^,K{x',x) h{x',t)) fail to describe 
the problem properly. As noted in the theoretical sec- 
tion one may then proceed systematically by studying 
the origins of the effective noise terms, e.g., on the basis 
of an assumption of Stokes' flow everywhere in the bulk. 
This is not only useful for the particular case of roughen- 
ing in imbibition but will also advance the understand- 
ing of the permeability properties, multi-phase flow and 
the role of capillary effects; for disordered porous media 
these are all very challenging problems (see among others 
j32> ..86..137,,.234. 292, 330J ). 

The fact that quantitative understanding of kinetic 
roughening should also involve considering the solu- 
tions to these problems again point to the complications 
present in imbibition phenomena. We might then con- 
clude that this merits a statistical physics approach, in 
which one eschews detailed considerations in favour of 
basic symmetry and scaling principles, but it is also im- 
portant to realise that imbibition represents an important 
meeting point for statistical physics and actual applica- 
tions to disordered media. 



B. Future issues for theorists 

The primary importance of imbition points to a num- 
ber of directions for theorists' future work. We now out- 
line briefly three separate topics, that stress quite differ- 
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ent aspects. 

Interface roughness: from the viewpoint of general 
interface dynamics it would of course be most interesting 
to develop models that correctly describe the experimen- 
tal results on interface roughness properties. One of the 
fundamental questions remaining open is of course the 
value (or possible values) of the roughness exponent x 
in a "most clear-cut" imbibition experiment. In the the- 
ory section the idea of "dynamical noise" was brought 
forward to explain some of the scaling regimes of the 
phase field-type models corresponding to constant fluid 
flow setups. This kind of reasoning is an analogy of the 
quenched-annealed cross-over in ordinary interface depin- 
ning (QEW or QKPZ), which results from the fact that 
the noise becomes decorrelated due to the finite velocity. 
This highlights the point that, as in local models, it is 
the noise and the symmetries of the effective interface 
dynamics that dictate the "critical" exponents - however 
with the proviso that the scaling range might be limited 
due to the intrinsic nature of imbibition. 

In this sense, attempts to quantitatively describe a 
variable-gap Hele-Shaw cell give rise to some optimism 
that the quenched features of the noise could be de- 
scribed, even quantitatively, in some simple cases. Lest 
it be forgotten we should emphasise that the theoret- 
ical and numerical results presented in Section IIIII are 
inspired by simple experiments and computational ease: 
no work exists so far on imbibition-like models in, say, 
three dimensions. In the context of various numerical 
(lattice Boltzmann or network) models this is a challeng- 
ing task. It is clear that the presence of an extra dimen- 
sion makes analytical work more challenging, and per- 
haps prese nts n ew scenarios as in the case of anisotropic 
depinning |32lj . It is also rather obvious that this would 
be of experimental and practical interest. 

Washburn behaviour: it is clear that in many in- 
stances of spontaneous imbibition Washburn-like scaling 
is not observed, and other (velocity) laws are obtained. 
This can ensue from the dynamic response of the medium 
(swelling) , from microscopic physics at the meniscus level 
(surfactant dynamics, inertial flow in transients), and 
from the fact that the liquid does not obey Stokes' flow 
conditions in general. Also, for e.g. very small capillary 
numbers it is possible that precursors (film flow) or the 
particular, peculiar pore structure play a role. We be- 
lieve that several of these effects merit separate studies, 
particularly in connection to their experimental effects on 
the roughness. The statistical physics viewpoint would 
be that such microscopic phenomena give rise to yet an- 
other length-scale: beyond this the effect of liquid conser- 
vation should dominate as usual. Now the question is to 
understand the dependence on time and on quantifiable 
measures of microscopic structure, of such scales. 

The coarse-grained description of these effects can also 
be improved; the phase field models described earlier do 
not do yet a very good job on the actual fluid dynamics 
level. Generalisations to account for this are of course 
possible. The scaling away from the micro-scale towards 



an effective large-scale description presents in the context 
of imbibition several challenges. At least three can be 
listed. In spontaneous dynamics, the time-scales change 
continuously as the interface slows down (recall that even 
without an asymptotic equilibrium dth ^ ^/h). Second, 
as noted the dynamics on the pore level can be very slow 
(coming from surfactant diffusion or induced dissipation 
at the contact line) or very fast (Bosanquet-flow). Third, 
the correlations in the pore structure and their size distri- 
bution can resist attempts to average, without even men- 
tioning time-dependent structural changes. In particular, 
also in relation to the actual kinetic roughening proper- 
ties, it is challenging to describe the dynamics of imbibi- 
tion when the interface advancement becomes slow, and 
stick-slip motion typical of general depinning/pinning 
transitions takes place. In the experimental Section and 
in the Introduction we pointed out evidences for devia- 
tions from Washburn-like approaches to a pinning height 
at late stages. These can in some cases arise from a 
change in the morphology of the interface from a com- 
pact (but rough) front, perhaps, and indeed involve fur- 
ther scales. Likewise, one should note that it should be 
evidently possible to construct phase-field -style theories 
that would allow to handle the smoothening and disper- 
sion of the front e.g. in the presence of initial wetting 
fluid saturations. 

Complete stochastic description: another promis- 
ing avenue for future work is to concentrate on issues 
other than those dealing with traditional "scaling expo- 
nents" . The dynamics of an interface in the presence 
of external (thermal) and quenched noise is of course 
a stochastic process and can be characterised in many 
ways. In the statistical physics community, much atten- 
tion has been recently paid to concep ts like persistence, 
or equivalently, first-return properties |l46lll65ll208ll20^ 
210, 229]. For example, the local interface height h{x,t) 
(after subtracting the average) and the local velocity of 
the interface v{x,t) are both stochastic variables. As in 
the simple cases of Brownian motion or ordinary random 
walks, it is then possible to explore the time-series of h, 
or e.g. return-to-origin statistics. In usual kinetic rough- 
ening models, which are most often in the steady-state, 
it follows that this defines a persistence exponent, sim- 
ply related to the exponent /3 of the two-point temporal 
correlation function if the dynamics is Markovian. The 
probability distributions or scaling functions of various 
statistical and fluctuating quantities, such as the velocity 
distribution P{v, L) (see Fig.l26|l or the roughness distri- 
bution in depinning problems of the QEW type P(w, F) 
|287j can are also similar alternatives descriptions. There 
has also been some recent interest in distributions such 
as "P(t')" as general signatures of non-equilibrium be- 
haviour and the possible universality therein [s^ . 

It has to be stressed however that the non-local dy- 
namics of imbibition does not necessarily maintain local 
scale-invariance, and thus many of these questions have 
to be re-evaluated. Also, one probably cannot expect an- 
alytical solutions as in the case of KPZ interfaces with 
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their deep connections to random matrices j267| . 

The scahng exponents of ordinary local interface equa- 
tions are also manifest in many other properties. Exam- 
ples are found in non-equilibrium steady states, as in the 
SOC ensemble 3], and in case the system is driven, e.g., 
in an oscillatory manner. In this case, one runs into con- 
cepts such as hysteresis, aging, and the general effect of 
"AC" driving fields 107. .25Q |. In the context of imbibi- 
tion it is clear that analogous scenarios are easy to set 
up in experiments by, e.g., varying the injection rate in 
a forced fluid flow test. 



C. Experimental and practical implications 

To finish the discussion on possible directions for fur- 
ther work we outline some suggestions for experimental- 
ists. Kinetic roughening and related questions should 
be studied in systems, where the micro-structure is un- 
der control — which is certainly not the case for the 
"classical" imbibition of water into paper. The often- 
used micro-channel networks are perhaps too idealised 
for such purposes. First of all, the average flow and its 
influence on the lateral length scale should be ad- 
dressed. Roughness exponents should be measured very 
carefully, not by a mere bona fide fit to various moments 
of height differences {\h{x + £,,t) - but also 

by analysing the structure factor and the probability dis- 
tribution functions of these quantities. 

The predictions of the phase field model can be applied 
directly also to the length-scales. For the <^x , Eq. H32() 
implies that a knowledge of the most important fluid pa- 
rameters and a characterisation of the porous medium 
can be used to estimate its range and effect s on scaling. 
Consider e.g. a Hele-Shaw cell, analogous to |l2l| . Using 
water and beads of size 0.4 mm, and varying the flow ve- 
locity between 10^® to 10^'^ m/s yields capillary numbers 
in the range 10~^ to 10^^. Assuming now that the typical 
pore size is of the same order of magnitude as the beads, 
this yields (with a surface tension of ^ 72.5 mJ/m ) ^x- 
values of about 1.5 mm for Ca = 10"^. In other words, 
a rather small length-scale if compared to typical kinetic 
roughening experiments. For ordinary paper the pores 
are smaller by only one order of magnitude. This shows 
that the roughness often seen in paper-based imbibition 
experiments may directly be due to the structural disor- 
der which is discrete on the sub-millimetre scales (recall 
that a typical fibre length is up to 2 to 4 mm) and does 
not follow any real "continuum description" . 

A central question for the dynamics is interface prop- 
agation by avalanches close to pinning. The complicated 
behaviour of simultaneous avalanches, in a system with 
non-local dynamics, is not understood well theoretically 
as already Fig. |^ above indicates. Possibly it may even- 
tually become feasible to explore the details of flu id dy- 
namical fiuctuations (pressure), as e.g. Ref. tSSlj hints. 
The question of pore-scale physics in this context and 
its relation to interface dynamics would seem to merit 



attention |l85j . 

Last we recall that quantitatively by far the largest 
amount of work has been done in detailed studies of phe- 
nomena relevant to all kinds of applications. Our under- 
standing is that this will remain so, since many of the 
phenomena listed in the experimental Section III can be 
examined more thoroughly: Pre-wetting layers along and 
inside pores, the influence of surfactants and solubles, 
partial saturation, the flow and transport in the bulk, 
and in particular behaviour with non-Newtonian liquids. 
Advances in various experimental techniques (X-ray to- 
mography for pore structure determination, NMR with 
high enough temporal resolution for dynamical measure- 
ments) should also provide with new results. As men- 
tioned in the Introduction, such methods are very close 
to getting to the level of imagining imbibition fronts with 
sufficient accuracy, and comparisons can then be made to 
the local pore structure including cross-correlations. 



D. Last outlook 

In this Review we have collected experimental and the- 
oretical studies of imbibition and tried to place them in 
a common framework of understanding the phenomenon. 
This can not be gained without including the lateral 
length-scale ^x appearing in the interface fluctuations, 
which is a central point in our theory Section. Also the 
scaling behaviour of interface fluctuations has to be stud- 
ied carefully. Different methods have to be compared 
for consistency, in particular most valuable information 
can be obtained from the structure factor S{k,t) and the 
probability distribution of height and velocity fiuctua- 
tions. 

To finish this work we mention a few directions that 
should be interesting to follow or which may very well 
become more important in the future than what one 
might conclude from our exposition of the field here. 
First of all, we hope that ideas about the nature of sys- 
tems where disorder and global conservation laws com- 
bine can benefit from the phase field lessons. One partic- 
ular example could be the description of vortex phases in 
supe rconductors and the coarse-grained theories thereof 

Second, we have not discussed the complications that 
ensue in the presence of more than two phases (three 
phase imbibition) except briefly in the context of surfac- 
tants. These systems are of obvious interest in some oil 
recovery scenarios, but have so far been analy sed mainly 
numerically in the literature [zE Hill Ell il82, 211i|. 
Here the pore-scale complications get even more mani- 
fold, and it is obvious that coarse-graining these into any 
kind of continuum theory is a challenge. Simi lar scenar - 
ios exist if non-Newtonian fluids are involved |l55l l327| : 
the rate-dependence of the liquid(s) will make itself man- 
ifest already inside single pores, and if one for instance 
considers the slow-down in spontaneous penetration, it is 
evident that the average flow of the wetting liquid is hard 
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to quantify. This area is however apt for much future de- 
velopment given the natural role of surfactants in many 
practical applications, and the fact that e.g. in petroleum 
industry -related scenarios even the basic const ituents - 
like crude - have very non-ideal properties |24l| . 

Another aspect not mentioned is the presence of si- 
multaneous transport processes (e.g. heat) during im- 
bibition (see J26J for the effect of thermal gradients). 
This can be complicated further, if some of the phases 
involved undergo chemical reactions (see [t^ and refer- 
ences therein) . The fingering properties of such processes 
have been investigated in Hele-Shaw cells, but imbibition 
-related conditions have not been considered in general. 
One may for instance consider a porous medium, where 
a more viscous fluid displaces another, while at the front 
chemical reactions take place. These may involve or even 
produce a surfactant, giving rise to a coupling to the con- 
tact angle. 

Third, we believe also that a sufficient understanding 
of imbibition will prove to be of interest in a variety of 
technological applications. Micro- and nano-machinery 
and chemically structured surfaces should provide ample 
examples of situations of practical and industrial inter- 



est |l05l 1275] , as imbibition into carbon nanotubes j318| . 
In such contexts the continuum description may fail due 
to the small scales that necessitate an atomistic treat- 
ment; consider for instance the physics of contact lines 
for a start. It should on the other hand be of interest to 
be able control imbibition properties (including dynam- 
ics) together with the permeability as hinted many times 
earlier, here. 
Acknowledgements : 

MJA would like to thank the Centre of Excellence pro- 
gram of the Academy of Finland for support and the 
SMC centre at La Sapienza, Rome, for hospitality. MD 
would like to thank the Canada Research Chair on 
Value Added Paper and Joe Aspler, Francois Drolet and 
Lyne Cormier for interesting discussions. MR was sup- 
ported by SFB 611 (Singuldre Phdnomene in mathema- 
tischen Modellen) of Deutsche Forschungsgemeinschaft 
and thanks Helsinki University of Technology for its hos- 
pitality. All authors wish to thank Ken R. Elder, Sami 
Majaniemi, and Tapio Ala-Nissila for collaboration on 
several works on imbibition theory. Also Janoinen Lohi 
is gratefully acknowedged for a long-standing and very 
fruitful collaboration. 



[I] R. Acevedo, A. Cuadrado, C. De la Torre, and S.M.D. de 
la Espina, Eur. J. Histochem. 46, 143 (2002) i 

[2] S. Akin, J.M. Schembre, S.K. Bhat, and A.R. Kovacek, J. 

Pet. Sci. Eng. 25, 149 (2000). 
[3] M.J. Alava, J. Phys. Cond. Mat. 14, 2353 (2002). 
[4] M. Alava and K. B. Lauritsen, Europhys. Lett. 53, 563 

(2001). 

[5] M.J. Alava, to appear in Advances in Condensed Mat- 
ter and Statistical Mechanics, Eds. E. Korutcheva and R. 
Cuerno, (Nova Science Publishers); cond-mat/0307668. 

[6] R. Albert, A.-L. Barabasi, N. Carle, and A. Dougherty, 
Phys. Rev. Lett. 81, 2926 (1998). 

[7] L.A.N. Amaral, A.-L. Barabasi, S.V. Buldyrev, S. Havlin, 
and H.E. Stanley, Phys. Rev. Lett. 72, 641 (1994). 

[8] D. Ambrosi and L. Preziosi, SIAM J. App. Math. 61, 22 
(2000). 

[9] D. Antonelh and A. Farina, Composites A30, 1367 (1999). 
[10] J.-Y. Arns,. C.H. Arns, A.P. Sheppard, R.M. Sok, M.A. 

Knackstedt, and W.V. Pinczewski, J. Pet. Sci. Eng. 39, 

247 (2003). 

[II] J. Asikainen, S. Majaniemi, M. Dube, and T. Ala-Nissila, 
Phys. Rev. E65, 052104 (2002). 

[12] J.S. Aspler, S. Davis, and M.B. Lyne, J. Pulp and Paper 

Sci. 13, J55 (1987). 
[13] J. Aspler, Nordic Pulp Paper Res. J. 1, 68 (1993). 
[14] J.C. Bacri, M. Rosen, and D. Salin, Europhys. Lett. 11, 

127 (1990). 

[15] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Lett. 59, 
381 (1987). 

[16] A.S. Balankin, A. Bravo-Ortega, and D. M. Matamoros, 

Phil. Mag. Lett. 80, 503 (2000). 
[17] A.S. Balankin, O. Susarrey, and J.M. Conzales, Phys. 

Rev. Lett. 90, 096101 (2003). 
[18] B.A. Baldwin and E.A. Spinier, J. Pet. Sci. Eng. 35, 23 



(2002). 

[19] A.-L. Barabasi, S.V. Buldyrev, S. Havhn, G. Huber, H.E. 
Stanley, and T. Vicsek, in Surface Disordering: Growth, 
Roughening and Phase Transitions, Eds. R. Jullien, J. 
Kertesz, P. Meakin and D.E. Wolf, Nova Science, Com- 
mack (1992). 

[20] A.-L. Barabasi and H.E. Stanley, Fractal Concepts in 
Surface Crowth, Cambridge University Press, Cambridge 
(1995). 

[21] H. Barkhausen, Physik. Zeitschr. 20, 401 (1919). 

[22] J. Bear and Y. Bachmat, Introduction to Modeling of 
Transport Phenomena in Porous Media, Kluwer Academic 
Publishers (1990). 

[23] A.Y. Beliaev and S.M. Hassanizadeh, Transp. Porous Me- 
dia 43, 487 (2001). 

[24] B. Berkowitz and R.P. Ewing, Surv. Geophys. 19, 23 
(1998). 

[25] M.G. Bernadiner, Transp. Porous Media 30, 251 (1998). 
[26] M.R. Bet, G. Goissis, S. Vargas, and H.S. Seliste-de- 

Araujo, Biomaterials 24, 131 (2003) 
[27] J. Bico and D. Quere, Europhys. Lett. 61, 348 (2003). 
[28] J. Bico, C. Tordeux, and D. Quere, Europhys. Lett. 55, 

214 (2001). 

[29] M. Blunt and P. King, Phys. Rev. A 42, 4780 (1990). 
[30] M. Blunt, M.J. King, and H. Scher, Phys. Rev. A 46, 
7680 (1992). 

[31] M.J. Blunt and H. Scher, Phys. Rev. E, 52 6387 (1995). 
[32] M.J. Blunt, Curr. Opin. Colloid Interface Sci. 6, 197 
(2001). 

[33] M.J. Blunt, M.D. Jackson, M. Piri, and P.H. Valvatne, 

Adv. Water Resour. 25, 1069 (2002). 
[34] J. Boiffin, L. Bruckler, and C. Aubry, Agronomic 3, 291- 

302 (1983). 

[35] D.T. Booth, J. Arid Environments 43, 91 (1999). 



49 



[36] C.H. Bosanquct, Phils. Mag. scr. 6 45, 525 (1923). 

[37] S.T. Bramwcll, P.C.W. Holdsworth, and J.F. Pinton, Na- 
ture 396, 552 (1998). 

[38] A.J. Bray, Adv. in Phys. 43, 357 (1994). 

[39] J. A. Bristow, Svensk Papperstidning 20, 645 (1971). 

[40] L. Bruckler, Agronomie 3, 213-222; ibid. 223-232 (1983). 

[41] R. Bruinsma and G. Aeppli, Phys. Rev. Lett. 52, 1547 
(1984). 

[42] S. Bryant, P.R. King, and D.W. Mellow, Transp. Porous 

Media 11, 53 (1993). 
[43] S.E. Buckley and M. Leverett, Trans. AIME 146, 107 

(1942). 

[44] Burgers, The nonlinear diffusion equation (1974). 
[45] S.V. Buldyrev, A.-L. Barabasi, F. Caserta, S. Havlin, 
H.E. Stanley and T. Vicsek, Phys. Rev. A 45 R8313 

(1992). 

[46] S.V. Buldyrev, A.-L. Barabasi, S. Havlin, J. Kertesz, 
H.E. Stanley, and H.S. Xenias, Physica 191 A, 220 (1992). 

[47] M. Cachile, R. Chertcoff, A. Calvo, M. Rosen, J.P. Hulin, 
and A.M. Cazabat, J. Colloid Interface Sci. 182, 483 

(1996). 

[48] J.W. Cahn and J.E. HilUard, J. Chem. Phys. 28, 258 
(1958). 

[49] J. Carmeliet, F. Descamps, and G. Houvenaghel, Transp. 
Porous Media 35, 67 (1999). 

[50] J.S. Ceballos-Ruano, T. Kupka, D.W. NicoU, J.W. Ben- 
son, M.A. loarmidis, C. Hansson, and M.M. Pintar, J. 
Appl. Phys. 91, 6588 (2002). 

[51] D.Y.C. Chan, B.D. Hughes, L. Paterson, and C. Sirakoff, 
Phys. Rev. A38, 4106 (1988). 

[52] D. Chang and M.A. loannidis, J. Coll. Int. Sci. 253, 159 
(2002). 

[53] P. Chauvc, P. Lc Doussal and K.J. Wicsc, Phys. Rev. 

Lett. 86, 1785 (2001). 
[54] K.S.A. Chen and L.E. Scriven, Tappi Journal 73, 151 

(1990). 

[55] Q. Chen, M.K. Gingras, and B.J. Balcom, J. Chem. Phys. 

119, 9609 (2003). 
[56] T. Chen, O. Chen, and P. Chen, J. Can. Pet. Technol. 

40, 29 2001 

[57] X. Chen, K.G. Kornov, Y.K. Kamath, and A.V. Neimark, 

Text. Res. Journ. 71, 862 (2001). 
[58] E. Chibowski, and F. Gonzalez-Caballero, Langmuir 9, 

330 (1993). 

[59] E. Chibowski, and L. Holysz, J. Adhcs. Sci. Technol. 11, 
1289 (1997). 

[60] T.W. Chiu, R.J. Wakeman, P.R. Harris, and O.F.J. 

Meuric, Chem. Engng. Res. & Design 74, 220 (1996). 
[61] M. Cil and J.C. Reis, J. Pet. Sci. Eng. 16, 61 (1996). 
[62] A. Clarke, T.D. Blake, K. Carruthers, and A. Woodward, 

Langmuir 18, 2080 (2002). 
[63] G.N. Constantinidcs and A.C. Payatakcs, Transp. Porous 

Media 38, 291 (2000) 
[64] R. Cuerno and M. Castro, Phys. Rev. Lett. 87, 236103 

(2001) . 

[65] H. Dahle and M. Celia Comput. Geosci. 3, 1 (1999). 
[66] E. Dana, and F. Skoczylas, Int. J. Multip. Flow 28, 1965 

(2002) . 

[67] H.T. Davis, R.A. Novy, L.E. Scriven, and P.G. Toledo, 

J. Phys. Cond. Mat. 2, SA457 (1990). 
[68] S.H. Davis and L.M. Hocking, Phys. Fluids 12, 1646 

2000. 

[69] P.G. De Gcnncs, In: Physics of Disordered Materials, 
eds. D. Adler et al. (Plenum, NY, 1985), p. 227. 



[70] P.G. dc Gcrmcs, Rev. Mod. Phys. 57, 827 (1985) 

[71] T. Dclkcr, D.B. Pcngra, and P.-z. Wong, Phys. Rev. Lett. 

76, 2902 (1996). 
[72] M. De Meijer, B. Van de Velde, and H. Militz, J. Coatings 

Tech. 914, 39 (2001). 
[73] A. De Wit, Phys. Fluids 15, 163 (2004). 
[74] M.M. Dias and A.C. Payatakes, J. Fluid Mech. 164, 305 

(1988). 

[75] R. Dickman, M. A. Mufioz, A. Vcspignani, and S. Zap- 
peri, Braz. J. Phys. 30, 27 (2000). 

[76] M.I.J, van Dijke and K.S. Sorbie, J. Pet. Sci. Eng. 39, 
201 (2003). 

[77] M.I.J, van Dijke and K.S. Sorbie, Paper SCA 2001-36, 
SCA Conference 2001, Edinburgh, September 2001. 

[78] F. Dommguoz, M.C. Gonzalez, and F.J. Cejudo, Planta 
215 727 (2002) 

[79] M.Z. Dong, F.A.L. Dullien, and J. Zhou, Transp. Porous 

Media 31, 213 (1998). 
[80] T. Dopier, A. Modaressi, and V. Michaud, Metall. Mat. 

Trans. B 31, 225 (2000). 
[81] A. Dougherty and N. Carle, Phys. Rev. E 58, 2889 

(1998). 

[82] M. Dube, M. Rost, K.R. Elder, M. Alava, S. Majaniemi, 

and T. AlarNissila, Phys. Rev. Lett. 83, 1628 (1999). 
[83] M. Dube, M. Rost and M. Alava, Eur. Phys. J. B 15, 691 

(2000). 

[84] M. Dube, M. Rost, K. R. Elder, M. Alava, S. Majaniemi 
and T. Ala-Nissila, Eur. Phys. J. B 15, 701 (2000). 

[85] M. Dube, S. Majaniemi, M. Rost, K. R. Elder, M. Alava, 
and T. Ala-Nissila, Phys. Rev. E 64, 051605 (2001). 

[86] F.A.L. Dullien, J. Porous Media 1, 29 (1998). 

[87] F.A.L. Dullien and M. Dong, J. Porous Media 5, 1 (2002). 

[88] K. Dcckolnick and G. Dziuk, Interfaces and Free Bound- 
aries 2, 341 (2000). 

[89] S.F. Edwards and D. Wilkinson, Proc. Roy. Soc. Loud. 
A (1982). 

[90] See, eg., D.E. Eklund, and P.J. Salminen, Tappi Journal 
69 (9), 116 (1986). 

[91] K. R. Elder, M. Grant, N. Provatas, and J. M. Kosterlitz, 
Phys. Rev. E 64, 21604 (2001). 

[92] A. Endruwcit, T. Luthy, and P. Ermanni, Polym. Com- 
pos. 23, 538 (2002). 

[93] B.J. Enquist, Plant, Cell & Environment 26, 151 (2003). 

[94] R.P. Ewing and S.C. Gupta, Water Resour. Res. 29, 3179 
(1993). 

[95] F. Family, K.C.B. Chan, J.G. Amar, in Surface Disor- 
dering : Growth, Roughening and Phase Transitions, R. 
JuUien, J. Kertesz, P. Meakin and D.E. Wolf, Nova Sci- 
ence, Commack (1992). 

[96] D.H. Fenwick and M.J. Blunt, Adv. Water Resour. 21, 
121 (1998). 

[97] J.M. Frey, P. Schmitz, J. Dufreche, and I.G. Pinheiro, 

Transp. Porous Media 37, 25 (1999). 
[98] R.E. Gagnon, G.D. Parish and D.W. Bousfield, Tappi 

Journal 84, 66 (2001). 
[99] V. Ganesan and H. Brenner, Phys. Rev. Lett. 81, 578 

(1998). 

[100] W.R. Carder and M.S. Mayhugh, Proc. Soil Sci. Soc. 

Amer.22, 187 (1958). 
[101] D. Geromichalos, F. Mugele, and S. Herminghaus, Phys. 

Rev. Lett. 89, 104503 (2002). 
[102] D. Geromichalos, personal communication (2003). 
[103] T. Giamarchi and S. Bhattacharya, In "High Magnetic 

Fields: Applications in Condensed Matter Physics and 



50 



Spectroscopy", p. 314, ed. C. Berthier et al., Springer- 
Verlag, 2002 (also cond-mat/01 11052). 
104] T. Gillespie, J. Coll. Science 13, 32 (1958), and ibid 14, 
123 (1959). 

105] N. Giordano and J.T. Cheng, J. Phys. Cond. Mat. 13, 
R271 (2001). 

106] M. Gladkikh and S. Bryant, Adv. Water Resour. 26, 

609 (2003). 

107] A. Glatz, T. Nattermann, and V. Pokrovsky Phys. Rev. 

Lett. 90, 047201 (2003). 
108] W. Gray and S. Hassanizadeh, Water Resour. Res. 27, 

1855 (1991). 

109] W.G. Gray and C.T. Miller, Phys. Rev. E61, 2150 

(2000). 

110] M.L.H. Gruwcl, B. Chatson, X.S. Yin, S. Abrams, Int. 

J. Food Sci. Tech. 36, 161 (2001). 
Ill] M.L.H. Gruwcl, X.S. Yin, M.J. Edney, S.W. Schroeder, 

A.W. MacGregor, S. Abrams, J. Agric. Food Chem. 50, 

667 (2002) 

112] J. D. Gunton, M. San Miguel, and P. Sahni, in Phase 
Transitions and Critical Phenomena, edited by C. Domb 
and J. L. Lebowitz (Academic Press, London, 1983), vol 
8, p. 267. 

113] B.S. Gupta, Tappi Journal 71, 147 (1988). 

114] L. L. Handy, Trans. Am. Inst. Mech. Eng. 219, 75 

(1960). 

115] B.I. Halperin and P.C. Hohenberg, Rev. Mod. Phys. 49, 
435 (1977). 

116] T. Halpin-Healy and Y.-C. Zhang, Phys. Rep. 254, 215 
(1995). 

117] A. Hamraoui and T. Nylander, J. Coll. Int. Sci. 250, 

415 (2002). 

118] A. Hamraoui, K. Thuresson, T. Nylander, and V. 

Yaminsky, J. Coll. Int. Sci. 226, 199 (2000). 
119] M. Hashemi, B. Dabir, and M. Sahimi, Aiche J. 45, 

1365 (1999). 

120] J. A. Hayashi and A. Soria, Aiche J. 47 1513, (2001). 
121] S. He, G.L.M.K.S. Kahanda, and P.-z. Wong, Phys. 

Rev. Lett. 69, 3731 (1992). 
122] A. Hornandcz-Machado, J. Soriano, A.M. Lacasta, M.A. 

Rodriguez, L. Ramirez-Piscina, and J. Ortm, Europhys. 

Lett. 55, 194 (2001). 
123] R. Hilfer, Phys. Rev. E 55, 5433 (1997). 
124] R. Hilfcr and H. Bcsscrcr, Physica B 279, 125 (2000). 
125] M. Hilpcrt, R. Glantz, and C.T. Miller, Transp. Porous 

Media 51, 267 (2003). 
126] L.M. Hirsch and A.H. Thompson, Phys. Rev. E50. 2069 

(1994) . 

127] K.T. Hodgson and J. C. Berg, J. Colloid Interface Sci. 

121, 22 (1988). 
128] V.K. Horvath, F. Family, and T. Vicsek, J. Phys. A 24, 

L25 (1991). 

129] V. K. Horvath, F. Family, and T. Vicsek, Phys. Rev. 

Lett. 67, 3207 (1991). 
130] V.K. Horvath and H.E. Stanley, Phys. Rev. E 52, 5166 

(1995) . 

131] J.Q. Hou, E.J. Kendall, G.M. Simpson, J. Exp. Bot. 48, 
683 (1997). 

132] R.W. Hoyland, Trans. BPBIF Symp. Fiber- Water In- 
teractions in Papermaking (Oxford), p. 557-579 (1976). 

133] Y.L. Hshieh and B.L. Yu, Text. Res. Journ. 62, 677 
(1992). 

134] Y.L. Hshieh, B.L. Yu, and M.M. Hartzell, Text. Res. 
Journ. 62, 697 (1992). 



135] G. Huber, M.H. Jensen, and K. Sneppen, Phys. Rev. 

E52, R2133 (1995). 
136] R.G. Hughes and M.J. Blunt, Transp. Porous Media 40, 

295 (2000). 

137] J.P. Hulin, Adv. Colloid Interface Sci. 49, 47 (1994). 
138] D.K. Jackson, M. Nicodemi, G. Perkins, N.A. Lindop, 

and H.J. Jensen, Europhys. Lett. 52, 210 (2000). 
139] P.I. Janscn and R.L. Ison, Seed Sci. Tech. 22, 435 (1994) 
140] D. Jasnow and J. Viiials, Phys. of Fluids 7, 747 (1996). 
141] G.R. Jerauld and S.J. Salter, Transp. Porous Media 5, 

103 (1990). 

142] H. Ji, M.O. Robbins, Phys. Rev. A 44, 2538 (1991). 
143] J.F. Joanny and L. Leger, Rep. Prog. Phys. 55, 431 

(1992). 

144] M. Jost and K.D. Usadel, Phys. Rev. B 54, 9314 (1996). 
145] M.A. Kader and S.C. Jutzi, J. Agronom. Crop Sci. 188, 
286 (2002) 

146] H. Kallabis and J. Krug, Europhys. Lett. 45, 20 (1999). 
147] M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 

56, 889 (1986). 
148] M. Kardar, Physics Reports 301, 85 (1998). 
149] K. Kawasaki and T. Ohta, Prog. Theo. Phys. 68, 129 

(1982). 

150] P. Kechagia, Y.C. Yortsos, and P. Lichtner, Phys. Rev. 

E64, 016315 (2001). 
151] D.A. Kessler, H. Levine, and Y.H. Tu, Phys. Rev. A43, 

4551 (1991). 
152] E. Kissa, Text. Res. Journ. 66, 660 (1996). 
153] M.A. Knackstcdt, A.P. Sheppard, and M. Sahimi, Adv. 

in Water Res. 24, 257 (2001). 
154] B. Koiller, H. Ji, M.O. Robbins, Phys. Rev. B 45, 7762 

(1992). 

155] L. Kondic, M. J. Shelley, and P. Palffy-Muhoray, Phys. 

Rev. Lett. 80, 1433 (1998). 
156] J. Koplik and T.J. Lasseter SPE J. 25, 89 (1985). 
157] J. Koplik and H. Levine, Phys. Rev. B32, 280 (1985). 
158] A. Koponen, D. Kandhai, E. Hellen, M. Alava, A. Hoek- 

stra, M. Kataja, K. Niskanen, P. Sloot, and J. Timonen, 

Phys. Rev. Lett. 80, 716 (1998). 
159] K.G. Korncv and A.V. Noirnark, J. Coll. Int. Sci. 235, 

101 (2001). 

160] J. Krug and P. Meakin, Phys. Rev. Lett. 66, 703 (1991). 
161] J. Krug and H. Spohn, in Solids Far From Equilibrium, 

Ed. C. Godreche, Cambridge University Press, Cambridge 

(1992). 

162] J. Krug, Phys. Rev. Lett. 72, 2907 (1994). 

163] J. Krug and L.-H. Tang, Phys. Rev. E50, 104 (1994). 

164] J. Krug, Adv. Phys. 46, 139 (1997). 

165] J. Krug, H. Kallabis, S.N. Majumdar, S.J. Cornell, A.J. 

Bray and C. Sire, Phys. Rev. E 56 2702 (1997). 
166] P.B.S. Kumar and D. Jana, Physica A 224, 199 (1996). 
167] M.C. Kuntz and J.P. Sethna, Phys. Rev. E 62 11699 

(2000). 

168] M. Kuntz, J. Van Mier, and P. Lavallee, Transp. Porous 

Media 43, 289 (2001). 
169] M. Kuntz and P. Lavallee, J. Phys D: Appl. Phys. 34, 

2547 (2001). 

170] T.H. Kwon, A.E. Hopkins, and S.E. O'Donnell, Phys. 

Rev. E 54, 685 (1996). 
171] L. Labajos-Broncano, M.L. Gonzalez-Martin, and J.M. 

Bruque, J. Adhes. Sci. Technol. 16, 1515 (2002). 
172] L. Labajos-Broncano, M.L. Gonzalez-Martin, J.M. 

Bruque, and CM. Gonzalez-Garcia, J. Coll. Int. Sci. 233, 

356 (2001). 



51 



173] L. Labajos-Broncano, M.L. Gonzalcz-Martm, J.M. 

Bmquc, CM. Gonzalcz-Garci'a, and B. Janczuk, J. Coll. 

Int. Sci. 219, 275 (1999). 
174] L. Labajos-Broncano, M.L. Gonzalez-Martin, B. 

Janczuk, J.M. Bruque, and C.M. Gonzalez- Garcia, J. Coll. 

Int. Sci. 211, 175 (1999). 
175] F. Lacombc, S. Zapperi and H. J. Herrmann, Phys. Rev. 

B63, 104104 (2001). 
176] M. Lago and M. Araujo, J. Coll. Int. Sci. 234, 35 (2001). 
177] M. Lago and M. Araujo, Physica A 289, 17 (2001). 
178] C.-H. Lam and V.K. Horvath, Phys. Rev. Lett. 85, 1238 

(2000). 

179] L.D. Landau and E.M. Lifshits, Gidrodinamika (Hydro- 
dynamics) 3rd. edition, Nauka, Moscow, (1986). 

180] J.S. Langer and L.A. Turski, Acta Metall. 25, 1113 
(1977). 

181] J.S. Langer, in Solids Far From Equilibrium, Ed. 
C. Godreche, Cambridge University Press, Cambridge 

(1992) . 

182] C. Larochc C, O. Vizika, and F. Kalaydjian, J. Pet. Sci. 
Eng. 24, 155 (1999). 

183] K. B. Lauritscn, P. Frojdh, and M. Howard, Phys. Rev. 
Lett. 81, 2104 (1998). 

184] P. Le Doussal, K.J. Wiese and P. Chauve, Phys. Rev. B 
66, 174201 (2002). 

185] R. Lenormand and C. Zarcone, Soc. of Petroleum En- 
gineers No. 13264, in: Proceedings of the 59th Ann. Tech. 
Conf. SPE (SPE, Richardson, TX, 1984). 

186] R. Lenormand, E. Touboul, and C. Zarcone, J. Fluid. 
Mech. 189, 165 (1988). 

187] R. Lenormand, J. Phys. Cond. Mat. 2, SA79 (1990). 

188] H. Loschhorn, Physica A 195, 324 (1993). 

189] H. Lcschhorn and L.-H. Tang, Phys. Rev. Lett. 70, 2973 

(1993) . 

190] H. Leschhorn and L.H. Tang, Phys. Rev. E 49, 1238 

(1994) . 

191] H. Leschhorn, Phys. Rev. E 54, 1313 (1996). 

192] H. Leschhorn, T. Nattermann, S. Stepanow, and L.-H. 
Tang, Ann. Physik 6, 1 (1997). 

193] S. Lotolicr, H.J. Lcuthcusscr, and Z. Rosas, J. Coll. Int. 
Sci. 72, 465 (1979). 

194] G. Leubner-Metzger, Planta 215, 959 (2002) effects on 
tobacco testa rupture and dormancy release 

195] A. Leventis, D.A. Verganclakis, M.R. Halsc, J.B. Web- 
ber and J.H. Strange, Transp. Porous Media 39, 143 
(2000) 

196] Y. Li, N.R. Morrow, and D. Ruth, J. Pet. Sci. Eng. 39, 
309 (2003). 

197] D.A. Lockington and J.Y. Parlange, J. Phys D: Appl. 

Phys. 36, 760 (2003). 
198] J.M. Lopez, M. A. Rodriguez, and R. Cuerno, Phys. Rev. 

E 56, 3993 (1997). 
199] J.M. Lopez, Phys. Rev. Lett. 83, 4594 (1999). 
200] P.J. Love, J.-B. Maillet, and P.V. Coveney, Phys. Rev. 

E 64, 061302 (2001). 
201] T.X. Lu, J.W. Biggar, and D.R. Nielsen, Water Resour. 

Res. 30, 3275 (1994). 
202] T.X. Lu, J.W. Biggar, and D.R. Nielsen, Water Resour. 

Res. 30, 3283 (1994). 
203] T.X. Lu, D.R. Nielsen, and J.W. Biggar, Water Resour. 

Res. 31, 11 (1995). 
204] P.E. Luner and D. Van Der Kamp, J. Pharmac. Sci. 90, 

348 (2001). 

205] B.B. Luokkala, S. Garoff, R.D. Tilton, and R.M. Suter, 



Langrnuir 17, 5917 (2001). 
[206] S.X. Ma, N.R. Morrow, and X.Y. Zhang, J. Pet. Sci. 

Eng. 18, 165 (1997). 
[207] J.-B. Maillet and P.V. Covcncy, Phys. Rev. E 62, 2898 

(2000). 

[208] S.N. Majumdar, Curr. Sci. (India), 77, 370 (1999). 
[209] S.N. Majumdar and A.J. Bray, Phys. Rev. Lett. 86, 
3700 (2001). 

[210] S.N. Majumdar and A.J. Bray, Phys. Rev. Lett. 91, 

030602 (2003). 
[211] V. Mani and K.K. Mohanty, SPE J. 3, 238 (1998). 
[212] V. Mani and K.K. Mohanty, J. Pet. Sci. Eng. 23, 173 

(1999). 

[213] J.A. Mann, L. Romero, R.R. Rye, and E.G. Yost, Phys. 

Rev. E 52, 3967 (1995). 
[214] P. Manneville and H. Chate, Physica D 96, 30 (1996). 
[215] A. Marmur and R.D. Cohen, J. Coll. Int. Sci. 189, 299 

(1997) . 

[216] J. Marro and R. Dickman Nonequilibrium Phase Tran- 
sitions in Lattice Models (Cambridge University Press, 
Cambridge, 1999). 

[217] J.A. Marsh, S. Garoff, and E.B. Dussan, Phys. Rev. 
Lett. 70, 2778 (1993). 

[218] M. MarsiU, A. Maritan, F. Toigo, and J.R. Banavar, 
Rev. Mod. Phys. 68, 963 (1996). 

[219] G. Martic, F. Centner, D. Seveno, D. Coulon, J. De 
Coninck, and T.D. Blake, Langrnuir 18, 7971 (2002). 

[220] N. Martys, M. Cieplak and M. O. Robbins, Phys. Rev. 
Lett. 69, 3193 (1991). 

[221] A. Maximenko and V.V. Kadet, J. Pet. Sci. Eng. 28, 
145 (2000). 

[222] M.B. McDonald, J. Sullivan, M.J. Lauer, Seed Sci. Tech. 
22, 79 (1994). 

[223] E. McEntyre, R. Ruan, R.G. Fulcher, Cereal Chem. 75, 
792 (1998) 

[224] P. Meakin, Fractals, scaling and growth far from equi- 
librium, Cambridge University Press, Cambridge, UK 

(1998) . 

[225] A. Medina, C. Pcrcz-Rosales, A. Pineda, F.J. Higuera, 

Rev. Mex. Fis. 47, 537 (2001). 
[226] A. Medina, A. Pineda, and C. Trevino, J. Phys. Soc. 

Jpn. 72, 979 (2003). 
[227] Y. Melean, D. Broseta, and R. Blossey, J. Pet. Sci. Eng. 

39, 327 (2003). 

[228] Y. Mclean, D. Broseta, A. Hasmy, and R. Blossey, Eu- 

rophys. Lett. 62, 505 (2003). 
[229] J. Merikoski, J. Maunuksela, M. Myllys, J. Timonen, 

and M. Alava, Phys. Rev. Lett. 90 024501 (2003) 
[230] V. Michaud and A. Mortensen, Compos. Pt. A-Appl. 

Sci. Manuf. 32, 981 (2001) 
[231] V. Michaud and J.A.E. Manson, J. Compos. Mater. 35, 

1150 (2001). 

[232] V. Michaud, R. Tornqvist, and J.A.E. Manson, J. Com- 
pos. Mater. 35, 1174 (2001). 

[233] M.-Carmen Miguel, J.S. Andrade Jr., and S. Zapperi, 
cond-mat/0304555; submitted to a special issue in the 
Brazilian Journal of Physics. 

[234] C.T. Miller, G. Christakos, P.T. Imhoff, J.F. McBridc, 
J.A. Pedit, and J.A. Trangcnstein, Adv. Water Resour. 
21, 77 (1998). 

[235] S.J. Mitchell, cond-mat/0210239. 

[236] I. Mitkov, D.M. Tartakovsky, and C.L. Winter, Phys. 

Rev. E58 R5245, (1998). 
[237] I. Mitkov, D.M. Tartakovsky, and C.L. Winter, Phys. 



52 



Rev. E61, 2152 (2000). 
[238] K. Mogcnscri and E.H. Stonby, Transp. Porous Media 
32, 299 (1998). 

[239] K.K. Mohanty, T. Davis, and L.E. Scriven, See. Pet. 

Eng. Reservoir Eng. 1, 113 (1987). 
[240] A. A. Moreira and J. S. Andrade Jr., Mendes Filho J. 

and S. Zapperi, Phys. Rev. B 66, 174507 (2002). 
[241] N.R. Morrow and G. Mason, Curr. Op. Coll. Int. Sci. 6, 

321 (2001). 

[242] T.E. Murnlcy, C.J. Radkc, and M.C. Williams, J. Col- 
loid Interface Sci. 109, 398 (1986), and ibid 109, 413 
(1986). 

[243] S. Mukherji and S.M. Bhattacharjee, Phys. Rev. Lett. 

79, 2502 (1997). 
[244] W.W. MuUins and R.F. Sekerka, J. Appl. Phys. 35, 444 

(1964). 

[245] M. Myllys, J. Maunuksela, M. Alava, T. Ala-Nissila, 
J. Merikoski and J. Timonen, Phys. Rev. E 64, 036101 
(2001). 

[246] T. Nagarninc and S. Miyazima, Fractals 4, 293 (1996). 
[247] O. Narayan and D. S. Fisher, Phys. Rev. B48, 7030 
(1993). 

[248] O. Narayan, Phys. Rev. E 62, R7563 (2000). 

[249] T. Nattermann, S. Stepanow, L.-H. Tang, and H. 

Leschhorn, J. Phys. (France) II 2, 1483 (1992). 
[250] T. Nattermann, V. Pokrovsky, and V.M. Vinokur Phys. 

Rev. Lett. 87, 197005 (2001). 
[251] C.J. Nederveen, Tappi Journal 77, 174 (1994). 
[252] M. Nicodemi and H.J. Jensen, J. Phys. A 34, Lll 

(2001). 

[253] K.J. Niskanen (ed.). Paper Physics, Fapet, Helsinki 
(1998). 

[254] K.J. Niskanen and M.J. Alava, Phys. Rev. Lett. 73, 
3475 (1994). 

[255] H.F. Nordhaug, M. CeUa, and H.K. Dahle, Adv. Water 

Resour. 26, 1061 (2003). 
[256] Z. Olami, I. Procaccia and R. Zeitak, Phys. Rev. E 49, 

1232 (1994). 

[257] P.E. Oren, S. Bakke, and O.J. Arntzen, SPE J. 3, 324 

(1998). 

[258] P.E. Orcn and S. Bakke, J. Pet. Sci. Eng. 39, 177 (2003). 
[259] X.R. Ouyang, T. van Voorthuysen, P.E. Toorop, 

H.W.M. Hilhorst, Int. Jour. Plant Sci. 163, 107 (2002). 
[260] M. Paczuski, S. Maslov and P. Bak, Phys. Rev. E 53, 

414 (1996). 

[261] E. Paune and J. Casademunt, Phys. Rev. Lett. 90, 
144504 (2003) 

[262] J.R.A. Pearson and P.M.J. Tardy, J. Non-Newt. Fluid 
Mech. 102, 447 (2002). 

[263] I. Pezron, G. Bourgaiii. and D. Quere, J. Colloid Inter- 
face Sci. 173, 319 (1995). 

[264] L.N. Pietrzak, J. Fregeau-Reid, B. Chatson, and B. 
Blackwell, Can. J. Plant Sci. 82, 512 (2002). 

[265] E.J. Pinthus and I.S. Saguy, J. Food Sci. 59, 804 (1994). 

[266] N. Poulin, P. Tanguy, J. Aspler, and L. Larrondo, Can. 
J. Chem. Eng. 75 949 (1997). 

[267] M. Prahofer and H. Spohn, Phys. Rev. Lett. 84, 4882 
(2000). 

[268] L. Preziosi, D.D. Joseph, and G.S. Beavers, Int. J. Mul- 
tiphase Flow 22, 1205 (1996). 

[269] N. Provatas, M.J. Alava, and T. Ala-Nissila, Phys. Rev. 
E54, R336 (1996). 

[270] D. Quere, Europhys. Lett. 39, 533 (1997). 

[271] D. Quere, E. Raphael, and J.Y. OUitrault, Langmuir 



15, 3679 (1999). 
[272] D. Rajagopalan, A. P. Aneja, and J.-M. Marchal, Text. 

Res. Journ. 71, 813 (2001). 
[273] E. Rame, J. Fluid Mech. 440, 205 (2001). 
[274] E.R. Rangel-German and A.R. Kovscek, J. Pet. Sci. 

Eng. 36, 45 (2002). 
[275] M. Rauscher and S. Dietrich, cond-mat/0303639 
[276] P. Reeves and M. Ceha, Water Resour. Res. 32, 2345 

(1996). 

[277] L.A. Richard, Physics (N.Y.) 1, 318 (1931). 
[278] E.K. Rideal, Philos. Mag. 44, 1152 (1922). 
[279] C.J. Ridgway and P.A.C. Gane, Colloid Surf. A 206, 
217 (2002). 

[280] C.J. Ridgway and P.A.C. Gane, Nordic Pulp Paper Res. 

J. 17, 119 (2002). 
[281] C.J. Ridgway, P.A.C. Gane, and J. Schoelkopf, J. Coll. 

Int. Sci. 252, 373 (2002). 
[282] R.F. Rodriguez, E. Sahnas-Rodriguez, J. A. Hayashi, A. 

Soria, and J.M. Zamora, Aiche J. 47, 1721 (2001). 
[283] S. Rods, J. Carmeliet, and H. Hens, Transp. Porous 

Media 52, 351 (2003). 
[284] J.G. Roof, Soc. Pet. Eng. J. 10, 85 (1970). 
[285] J. Rosinski, Am. Ink Maker 71, 40 (1993). 
[286] A. Rosso and W. Krauth, Phys. Rev. Lett. 87 (2001) 

187002. 

[287] A. Rosso, W. Krauth, P. Le Doussal, J. Vannimenus 
and K.J. Wiese, cond-mat/0301464 (2003). 

[288] M.A. Rubio, C. Edwards, A. Dougherty, and J. P. Gol- 
lub, Phys. Rev. Lett. 63, 1685 (1989). 

[289] R.R. Rye, J.A. Mann, and E.G. Yost, Langmuir 12, 555 
(1996). 

[290] M. Sacande, E.A. Golovina, A.C. van Aelst, F.A. Hoek- 

stra, J. Exp. Bot. 52, 919 (2001). 
[291] P.G. Saffman and G.I. Taylor, Proc. Roy. Soc. London 

A 245, 312 (1958). 
[292] M. Sahimi, Rev. Mod. Phys. 65, 1393 (1993). 
[293] E. J. Samuelsen, O.W. Gregersen, P.J. Houen, T. Helle, 

C. Raven, and A. Snigirev, J. Pulp Pap. Sci. 27, 50. 

(2001). 

[294] K. Saripalli, H. Kim, P. Rao, and M. Annablc, Environ. 

Sci. Technol. 31, 932 (1997). 
[295] A. E. Scheidegger, The Physics of Flow thought Porous 

Media, 3rd Ed., University of Toronto Press, Toronto 

(1974). 

[296] J.M. Schcmbrc and A.R. Kovscek, J. Pet. Sci. Eng. 39, 
159 (2003). 

[297] J. Schmittbuhl, J. Vilotte, and S. Roux, Phys. Rev. 

E51, 131 (1995). 
[298] J. Schoelkopf, P.A.C. Gane, C.J. Ridgway, and G.P. 

Matthews, Colloid Surf. A 206, 445 (2002) 
[299] J. Schoelkopf, P.A.C. Gane, C.J. Ridgway, and G.P. 

Matthews, Nordic Pulp Paper Res. J. 15, 422 (2000). 
[300] J. Schoelkopf, C.J. Ridgway, P.A.C. Gane, G.P. 

Matthews, and D.C.J. Spielmann, J. Coll. Int. Sci. 227, 

119 (2000). 

[301] D.R. Schuchardt and J.C. Berg, Wood Fib. Sci. 23, 342 

(1991). 

[302] T.J. Senden, M.A. Knackstedt, and M.B. Lyne, Nordic 
Pulp Paper Res. J. 15, 554 (2000). 

[303] E.T. Scppala, V. Petaja, and M.J. Alava, Phys. Rev. 
E58, R5217 (1998). 

[304] A. Siebold, M. Nardin, J. Schultz, and A. WalUser, Col- 
loid Surf. A 161, 81 (2000) 

[305] D.E. Smiles, Chem. Eng. Sci. 53, 2211 (1998). 



53 



[306] K. Snoppon, Phys. Rev. Lett. 69, 3539 (1992). 
[307] K. Sncppcn and M.H. Jensen, Phys. Rev. Lett. 71, 101 
(1993). 

[308] R.M. Sok, M.A. Knackstedt, A.P. Sheppard, W.V. 

Pinczewski, W.B. Lindquist, A. Venkatarangan, and L. 

Paterson, Transp. Porous Media 46, 345 (2002). 
[309] J.L. Sommer and A. Mortensen, J. Fluid. Mech. 311, 

193 (1996). 

[310] K.S. Sorbie, Y.Z. Wu, and S.R. McDougall, J. Coll. Int. 

Sci. 174, 289 (1995). 
[311] J. Soriano, J.J. Ramasco, M.A. Rodriguez, A. 

Hernandez-Machado and J. Ortm, Phys. Rev. Lett. 89, 

026102 (2002). 

[312] J. Soriano, J. Ortm and A. Hernandez-Machado, Phys. 

Rev. E 66 031603 (2002). 
[313] J. Soriano, J. Ortm and A. Hernandez-Machado, cond- 

mat/0208432. 

[314] D. Sornette, Critical Phenomena in Natural Sciences 

(Springer, Berlin, 2000). 
[315] K. Stoev, E. Rame, T. Leonhardt, and S. Garoff, Phys. 

Fluids 10, 1793 (1998). 
[316] K. Stoev, E. Rame, and S. Garoff, Phys. Fluids 11, 3209 

(1999). 

[317] J.P. Stokes. D.A. Weitz, J.P. Gollub, A. Dougherty. 

M.O. Robbins, P.M. Chaikin, and H.M. Lindsay, Phys. 

Rev. Lett. 57, 1718 (1986). 
[318] S. Supple and N. Quirke, Phys. Rev. Lett. 90, 214501 

(2003). 

[319] J. Szekely, A.W. Neumann, and Y.K. Chuang, J. Colloid 

. Interface Sci. 35, 273 (1971). 
[320] A. Takahashi, M. Haggkvist, and T.Q. Li, Phys. Rev. E 

56, 2035 (1997). 
[321] L.-H.Tang, M. Kardar,and D. Dhar, Phys. Rev. Lett. 
74, 920 (1995). 

[322] A. Tanguy, M. Gounelle, and S. Roux, Phys. Rev. E 58, 

1577 (1999). 
[323] K. E. Thompson, Aiche J. 48 1369, 2002. 
[324] C. Thorenz, G. Kosakowski, O. Kolditz, and B. 

Berkowitz, Water Resour. Res. 38, 1069 (2002). 
[325] F. Tiberg, B. Zhmud, K. Hallstensson, and M. von Bahr, 

Phys. Chem. Chem. Phys. 2, 5189 (2000). 
[326] CD. Tsakiroglou and A.C. Payatakes, Adv. Water Re- 
sour. 23, 773 (2000). 
[327] CD. Tsakiroglou, M. Theodoropoulou, V. Karoutsos, 

D. Panauicolaou, and V. Sygouni, J. Colloid Interface Sci. 

267, 217 (2003). 
[328] G. C. Tzimas, T. Matsuura, D.G. Avraam, W. Van der 

Brugghen, G.N. Constantinides, and A.C. Payatakes, J. 



Colloid Interface Sci. 189, 27 (1997). 
[329] G.C. Tzimas, T. Matsuura, D.G. Avraam, W. Van Der 

Brugghen, G.N. Constantinides, A.C. Payatakes, J. Coll. 

Int. Sci. 189, 27 (1997). 
[330] O. van Genabeek and D.H. Rothman, Annu. Rev. Earth 

Planet. Sci. 24, 63 (1996). 
[331] S. C. van der Marck and J. Glas, Eur. J. Mech. B-Fluids 

16, 681 (1997). 
[332] W. van Saarloos, Phys. Rep. 301, 9 (1998). 
[333] A. Vernhet, M.N. Bellon Fontaine, J.M. Brillouet, E. 

Roesink, and M. Moutounet, J. Membr. Sci. 128, 163 

(1997). 

[334] A.M. Vidales, J.L. Riccardo, and G. ZgrabUch, J. Phys. 

D Appl. Phys. 31, 2861 (1998). 
[335] O. Vizika, D.G. Avraam, and A.C. Payatakes, J. Colloid 

Interface Sci. 165, 386 (1994). 
[336] H.J. Vogel and K. Roth, J. Hydrology 272, 95 (2003) 
[337] M.E.P. Walinder and D.J. Gardner, J. Adh. Sci. Tech. 

13, 1363 (1999). 
[338] N.C Wardlaw and Li Yu Transp. Porous Media 3, 17 

(1988) 

[339] E.W. Washburn, Phys. Rev. 17, 273 (1921). 
[340] G.J. Weir and S.P. White, Transp. Porous Media 25, 79 
(1996). 

[341] D. Wildenschild, J.W. Hopmans, C.M.P. Vaz, M.L. 
Rivers, D. Rikard, and B.S.B. Christensen, J. Hydrology 
267, 285 (2002) 

[342] D. Wilkinson, Phys. Rev. A 34, 1380 (1986). 

[343] Y.S. Wu, K. Pruess, and P.A. Witherspoon, Transp. 
Porous Media 6, 115 (1991). 

[344] E.F. Wuebker, R.E. Mullen, and K. Koehler, Crop Sci- 
ence 41, 1857 (2001). 

[345] S. Zapperi, A. A. Moreira and J. S. Andrade Jr., Phys. 
Rev. Lett. 86, 3622 (2001). 

[346] S. Zapperi, J. S. Andrade Jr., and A. A. Moreira, cond- 
mat/0311034. 

[347] X.Y. Zhang, N.R. Morrow, and S.X. Ma, SPE Reserv. 

Engng. 11, 280 (1996) 
[348] D.X. Zhang, R.Y. Zhang, S.Y. Chen, and V.E. Soil, 

Geophys. Res. Lett. 27, 1195 (2000). 
[349] D. Zhou, L. Jia, J. Kamath, and A.R. Kovscek, J. Pet. 

Sci. Eng. 33, 61 (2002). 
[350] O. Zik, E. Moses, Z. Olami, and I. Webman, Europhys. 

Lett. 38, 509 (1997). 
[351] O. Zik, T. Kustanovicli, E. Moses, and Z. Olami, Phys. 

Rev. E 58, 689 (1998). 



